close

Вход

Забыли?

вход по аккаунту

?

A numerical study of microwave detection of a malignant tissue

код для вставкиСкачать
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A NUMERICAL STUDY OF MICROWAVE
DETECTION OF A MALIGNANT
TISSUE
The members o f the Committee approve the master’s
thesis of Luis Manuel Camacho-Velazquez
Saibun Tjuatja
Supervising Professor
Jonathan Bredow
Hanli Liu
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Copyright © by Luis Manuel Camacho-Velazquez 2004
All Rights Reserved
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
A NUMERICAL STUDY OF MICROWAVE
DETECTION OF A MALIGNANT
TISSUE
by
LUIS MANUEL CAMACHO-YELAZQUEZ
Presented to the Faculty of the Graduate School of
The University of Texas at Arlington in Partial Fulfillment
of the Requirements
for the Degree of
MASTER OF SCIENCE IN ELECTRICAL ENGINEERING
THE UNIVERSITY OF TEXAS AT ARLINGTON
May 2004
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
UMI N um ber: 1421247
INFORMATION TO USERS
The quality of this reproduction is dependent upon the quality of the copy
submitted. Broken or indistinct print, colored or poor quality illustrations and
photographs, print bleed-through, substandard margins, and im proper
alignm ent can adversely affect reproduction.
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if unauthorized
copyright material had to be removed, a note will indicate the deletion.
UMI
UMI Microform 1421247
Copyright 2004 by ProQuest Information and Learning Company.
All rights reserved. This microform edition is protected against
unauthorized copying under Title 17, United States Code.
ProQuest Information and Learning Company
300 North Zeeb Road
P.O. Box 1346
Ann Arbor, Ml 48106-1346
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ACKNOmEDGEMENTS
I would like to express my sincere appreciation to my supervising professor, Dr.
Saibun Tjuatja, for Ms guidance, assistance and stimulating discussion throughout my
Master of Science degree studies. My success in tMs M.S. program is largely due to his
time, encouragement, and practical insights. I would also like to thank the members of
the thesis committee, Dr. Jonathan Bredow, and Dr. Hanli Liu for their advices, as well
as for allocating their precious time in participating in my thesis defense examination.
As a faculty of Universidad Autonoma de Nuevo Leon, Mexico, I would also
like to thank all of my colleagues for their support and encouragement. Especially to
Jose A. Gonzalez, Jose L. Castillo, Castulo Vela, Rogelio Garza, and Eugenio Lopez.
Furthermore, I would like to share the satisfaction of this degree with my
classmates and excellent friends, Jose Mireles Jr., Humberto Ochoa, Tse-Min Chen, and
Cesar Heyaime for their discussions, help and generosity.
Finally, I would like to express my love and appreciation to my wife, Fina, and
my daughters, Estefania and Ximena, for their patience, understanding and M l support.
I also thank my parents, Jose Delfino Camacho, from whom I find out that learning is a
lifetime process, and Luvia Velazquez de Camacho, for her prayers and love. Last, but
not least, I want to thank my family and family-in-law, for their support. I dedicate this
work to all of them.
November 24, 2003
iv
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ABSTRACT
A NUMERICAL STUDY OF MICROWAVE
DETECTION OF A MALIGNANT
TISSUE
Publication No.
Luis Manuel Camacho-Velazquez, M.S.
The University of Texas at Arlington, 2004
Supervising Professor: Saibun Tjuatja
The physical basis for tumor detection with microwave imaging is the contrast
in the dielectric properties of background and malignant tissues.
In this research a
numerical study using the Finite-Difference Time-Domain method has been carried out
to analyze the use of microwave for lung cancer detection.
A computational model of a human torso slice was generated, and simulations
were conducted using several frequencies for several sizes and locations of a lung
cancerous tissue.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The tests were conducted by simulating a plane wave impinging on the torso
slice, and the analysis of the scattered field at different frequencies showed that 3 GHz
appears to be the optimum frequency for lung cancer detection.
A time series analysis is also presented as an approach to find the optimum
time, during the incidence of a plane wave, at which the presence of the tumor is more
discernible.
VI
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS...............................................
iv
ABSTRACT ............
v
LIST OF ILLUSTRATIONS
.....
ix
LIST OF TABLES .................
xv
Chapter
1. INTRODUCTION
.....
1
2. LUNG CANCER INFORMATION
................
6
3. FDTD METHOD DESCRIPTION.........................................................
17
4. COMPARISON OF FDTD RESULTS WITH THEORETICAL
AND PUBLISHED RESULTS...
..............
50
5. GENERATION OF THE COMPUTATIONAL MODEL
62
6. SIMULATIONS, RESULTS AND ANALYSES.......
7. CONCLUSIONS AND FUTURE W ORK
..........
.....
76
.....
106
Appendix
A. APPROXIMATION FOR DERIVATIVES
...................
B. GRAPHS OF DIELECTRIC PROPERTIES OF BIOLOGICAL
TISSUES FOR FREQUENCY RANGE 1 - 30 G H Z ........ .
C. MATLAB PROGRAM “TORSO_DIMENSIONS_IN_CELLS”
109
112
...... .
116
D. MATLAB PROGRAM “IMAGE2C0NSTANTS” .....................................
123
E. FDTD SIMULATION CODE
127
.....
vii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
REFERENCES.........................................................
146
BIOGRAPHICAL INFORMATION.........................................................................
152
V lll
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
LIST OF ILLUSTRATIONS
Figure
Page
2.1
Respiratory system..........................................................................................
7
2.2
a) Normal cells; and b) cells forming a tumor
8
2.3
CT scan showing a tumor........
3.1
Positions of field components. The E-components are in
the middle of edges and the H-components are in the
.....
center of the faces.
18
3.2
Leapfrog time integration.
22
3.3
Plane-wave propagation along the diagonal of a
two-dimensional square-cell ( A x x A x ) mesh is equivalent
to one-dimensional propagation along a mesh with cell size
...............................
......
12
......
Ax/V 2 , thus accounting for the two-dimensional stability
limit of cAt < k x l s j l
.....
25
3.4
The FDTD unit cell for the two-dimensional TM case.
................................. 34
3.5
The FDTD unit cell for the two-dimensional TE case.
.................................. 35
3.6
TOta! field/scattered field of the two-dimensional
problem space. ................................................................................................ 37
3.7
Total field/scattered field of the two-dimensional
problem space.
3.8
3.9
.........
38
2-D simulation of a plane wave propagating in free space.
The pulse is generated at one end and subtracted out at
the other end.
.....
39
Results of a simulation to show qualitatively the effectiveness
of the absorbing boundary condition. .....
44
IX
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.10
Phase speed (relative to c) versus number of mesh points
per wavelength for three different stability factors. ....................................... 46
3.11
Group speed (relative to c) versus number of mesh points per
wavelength for three different stability factors............................................... 47
3.12
Variation of the numerical phase velocity with wave
propagation angle in a two-dimensional FDTD grid
for three grid resolutions. ............................................................................... 49
4.1
Configuration of FDTD space to verify the reflection
coefficient for an EM plane wave traveling from free
space to a dielectric medium with S r- 4 and a = 0, ....................................... 51
4.2
EM pulse traveling from free space to a dielectric
medium with
= 4. As expected, the amplitude of
the transmitted wave is 0.66, while the reflected
wave is 0.33. .............
52
EM sinusoidal wave traveling from free space to a
dielectric medium with s ^ - 4 ,
and a = QS I m.
The wavelength inside the medium changes accordingly
to the permittivity.
.....
54
Configuration of the simulation o f a plane wave impinging
on the dielectric cylinder.
.................
55
4.3
4.4
4.5
4.6
Uniform plane wave incident on a dielectric circular
cylinder.
......
Simulation of a uniform plane wave incident on a dielectric
.....
circular cylinder.
56
57
4.7
Comparison of the FDTD results vs. Bessel function
expansion results along the propagation center axis of
a cylinder at three frequencies. The cylinder is 20 cm
in diameter and has a dielectric constant of 30 and a
conductivity of 0.3 S/m. ................................................................................. 58
4.8
Model and characteristics of the tissues in a human head.
4.9
Electric field distribution |Ez| in the biological structure
according to (a) Finite-Element Procedure and (b) FDTD. ............................ 61
.............
X
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
60
5.1
5.2
5.3
5.4
Discretization of a surface. For illustrative purposes the
cell size was not small enough to make staircasing
problem more e v id e n t
...... ............................. .......................... .
63
Image of a cross section of thorax from Stanford
University Medical Media and Information
Technologies Web Page
.....
64
Simplified version of a thorax cross section. These
organs within the thorax are: 1) lungs; 2) blood, which
fills the atria and the ventricles; 3) bones (ribs and
backbone); 4) heart wall; 5) interstitial fluid; 6)
esophagus; 7) intercostal muscles; and 8) thoracic wall,
.......
which includes the skin, the fat, and the m uscles.
65
Flow diagram of generation of the computational
torso model.
......
66
......
67
5.5
Dimensions o f the simulated thorax slice.
5.6
Lungs dielectric properties for the frequency range 1-30 GHz.
Values taken from the Institute for Applied Physics, Italy;
based on parametric model developed by C. Gabriel et al.
.........
68
5.7
Tumor dielectric properties for the frequency range 1-30 GHz. .................... 70
5.8
Values of relative permittivity ( s ^ ) for each of the tissues
for frequency = 5 G H z.
.....
73
........
74
5.9
Array containing the identifiers for a 5-GHz analysis.
6.1
Dielectric properties of lung and tumor. ........................................................ 77
6.2
Diagram of the simulation of an EM wave impinging
on the torso model. ......................................................................................... 78
6.3
Positions of the tumor for the study a) Position 1, at the
outer edge o f the right lung, and b) Position 2, at the inner
edge of the left lung. .....
6.4
Amplitude of the electrical field Ez read at the screen for
the torso model without tumor. The simulated frequencies
xi
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
79
are 1, 5,10,14, 18,and 28GHz.
6.5
.......
Amplitude of the electrical field Ez read at the screen for
the torso mode! without tumor. The simulated frequencies
are 1, 2, 3,4, 5, and 6GHz.
......
80
81
6.6
Diagram of case 1. EM plane wave impinging on the torso
slice with a 1-cm diameter tumor located at outer edge
of right lung. ................................................................................................... 83
6.7
Difference signal for case 1. EM plane wave impinging on
the torso slice with a I-cm diameter tumor located at outer
edge of right lung.
.....
84
Diagram of case 2. EM plane wave impinging on the torso
slice with a 3-cm diameter tumor located at outer edge
of right lung.
.....
85
Difference signal for case 2. EM plane wave impinging on
the torso slice with a 3-cm diameter tumor located at outer
edge of right lung...............................................................................
86
Diagram of case 3. EM plane wave impinging on the torso
slice with a 5-cm diameter tumor located at outer edge
of right lung
.....
87
Difference signal for case 3. EM plane wave impinging on
the torso slice with a 5-cm diameter tumor located at an
outer edge of right lung.
.............
88
6.8
6.9
6.10
6.11
6.12
Diagram of case 4. EM plane wave impinging on the torso
slice with a 1-cm diameter tumor located at inner edge
of left lung. ..................................................................................................... 89
6.13
Difference signal for case 4. EM plane wave impinging on
the torso slice with a 1-cm diameter tumor located at inner
edge of left lu n g . ........
90
6.14
Diagram of case 5. EM plane wave impinging on the torso
slice with a 3-cm diameter tumor located at inner edge
of left lung. ..................................................................................................... 92
6.15
Difference signal for case 5. EM plane wave impinging on
the torso slice with a 3-cm diameter tumor located at inner
xii
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
edge o f left lu n g .
6.16
6.17
......
93
Diagram of case 6. EM plane wave impinging on the torso
slice with a 5-cm diameter tumor located at inner edge
of left limg.
.....
94
Difference signal for case 6. EM plane wave impinging on
the torso slice with a 5-cm diameter tumor located at inner
edge of left lung.
......................
95
6.18
Difference signal (AEz) for 3-GHz steady state analysis.
The tumor is located on outer edge of right lung. .......................................... 96
6.19
Difference signal (AEz) for 3-GHz steady state analysis.
The tumor is located on inner edge of left lung.
6.20
.....
Diagram of the time series analysis. The values at the
screen are shown at every time step
96
98
6.21
Time series of scattered AEz at the screen for a 1-GHz
incident wave. A 1-cm tumor is located at outer edge
of right lung.......... ............................................................................... ..............98
6.22
Time series of scattered AEz at the screen for a 1-GHz
incident wave. A 5-cm tumor is located at outer edge
of right lung.
..... .................................................................................... 99
6.23
Time series of scattered AEz at the screen for a 1-GHz
incident wave. A 1-cm tumor is located at inner edge
of left lung......................
100
Time series of scattered AEz at the screen for a 1-GHz
incident wave. A 5-cm tumor is located at inner edge
of left lung.
.....
101
6.24
6.25
Time series of scattered AEz at the screen for a 3-GHz
incident wave. A 1-cm tumor is located at outer edge
of right lung. ................................................................................................... 101
6.26
Time series of scattered AEz at the screen for a 3-GHz
incident wave. A 1-cm tumor is located at outer edge
of right lung. ................................................................................................... 102
xni
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
6.27
6.28
6.29
6.30
Time series of scattered AEz at the screen for a 3-GHz
incident wave. A 1-cm tumor is located at inner edge
of left lung.
.....
102
Time series of scattered AEz at the screen for a 4-GHz
incident wave. A 1-cm tumor is located at outer edge
of right lung.
......
103
Time series of scattered AEz at the screen for a 4-GHz
incident wave. A 1-cm tumor is located at inner edge
of left lung.
.........
104
Time series of scattered AEz at the screen for a 6-GHz
incident wave. A 1-cm tumor is located at outer edge
o f right lung.
.....
104
XIV
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
LIST OF TABLES
Table
Page
2.1 Characteristics of Small Cell Lung Cancer..................................................... 10
2.2
Characteristics of Squamous Cell Lung Cancer............................................
2.3
Characteristics of Adenocarcinoma........
10
.......................
11
2.4 Characteristics of Large Ceil Carcinoma........................................................ 12
3.1
Location o f the different fields in the Yee cell.................................
19
5.1 Values of relative permittivity and conductivity for the
tissues at 1 GHz............................................................
6.1
Description of cases to analyze
7
....................................
6.2 Valuesof maximum peak of the difference signal for case 1
82
..............
6.3 Valuesof maximum peak of the difference signal for case 2 ......................... 86
6.4 Values of maximum peak of the difference signal for case3 ......................... 88
6.5 Values of maximum peak of the difference signal for case4......................... 91
6.6 Values of maximum peak of the difference signal for case 5 ......................... 93
6.7 Valuesof maximum peak of the difference signal for case 6......................... 95
XV
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
84
CHAPTER 1
INTRODUCTION
Lung cancer is now the most common form of cancer diagnosed in the United
States and a major cause of death. Lung cancer accounts for 14% of all cancers and
28% of all cancer deaths [!]. Besides, lung cancer is known as one of the most difficult
cancer to cure, and the number of deaths that it causes is usually increasing, as the
disease is detected at its evaluate stage. Because of the large size of the lungs, cancer
may grow for many years, undetected, without causing suspicion. Lung cancer can
even spread outside the lungs without causing severe symptoms.
About 75% of lung cancer patients experience a cough prior to diagnosis. Other
symptoms that are often ignored or explained away as signs of aging are dyspnea
(difficult breathing) and fatigue [2]. These are all reasons why lung cancer often goes
undetected until it is quite advanced.
Early diagnosis of lung cancer is not common because screening for lung cancer
is very rarely done. Even though chest x-rays can detect small tumors and CT scans can
detect even smaller tumors, these tests are generally requested only if symptoms are
present.
The tremendous toll that cancer takes, combined with the risks associated with
the limited sensitivity and specificity of conventional techniques, motivate the search
1
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
for alternative teclmologies for early cancer detection. One such alternative is
microwave imaging.
Microwave imaging for medical applications has been of interest for many years
[3]. Microwave imaging is defined as “seeing” the internal strachire of an object by
means of electromagnetic fields at microwave firequencies (300 MHz - 30 GHz) [4],
Microwave images are maps of the electrical property distributions in the body.
The electrical properties of various tissues may be related to their physiological
situation. For example, the properties of tissues change with temperature. Monitoring
hyperthermia has been proposed as one application of microwave imaging [5].
Hyperthermia is the application of heat to tissue. In this case, the changing electrical
properties indicate the successful deposition of heat in the tissue of interest.
Disease may cause other changes in electrical properties.
There is some
evidence of change in the properties of cancerous tissues when compared to normal
tissues.
Hence, the motivation for developing a microwave imaging technique for
detecting cancer is the significant contrast in the dielectric properties at microwave
frequencies of normal and malignant tissues suggested by published measured data [6].
Cancer detection with microwave imaging is based in this contrast in electrical
properties.
Recently, microwave imaging for breast cancer detection has interested
many researchers [7-11].
Concerning the lung cancer, most of the conducted research to detect a
cancerous tissue in the lung area is performed by processing either X-ray images or MR
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
images [12],
Nevertheless, with imaging tedmiques such as x-ray computed
tomography (CT) or magnetic resonance imaging (MEJ), it is not always possible to
distinguish malignant from norma! tissue [13]. To the best of our knowledge there is
not a significant research using microwaves to detect lung cancer.
Several advantages are expected from microwave imaging. First, it is totally
noninvasive. Second, the energy of “photons” in the microwave region is small enough
to avoid ionization effects, which are found with X-ray tomography. Microwave cancer
detection systems are expected to operate with power levels one or two orders of
magnitude below those of a cellular telephone. Evaluating SAR for specific systems
with computer models ensures safety and compliance with standards. Consequently,
microwave imaging is not expected to represent a health risk to the patient.
Microwave cancer detection is also expected to be less expensive than methods
such as MRI and nuclear medicine.
This is because microwave equipment costs a
fraction o f the equipment needed for MRI and nuclear medicines. The imaging process
is also anticipated to be very rapid, sensitive and specific.
The key to sensitivity,
specificity and the ability to detect small tumors is the anticipated electrical property
contrast between malignancies and normal tissues.
Specifying whether a suspicious
area is malignant or benign may be possible with a contrast in the electrical properties
of these tissues.
In order to simulate the study presented in this thesis a numerical technique was
applied. Numerical techniques are used to perform a simulation of electromagnetic
wave propagation in different media [14]. For this purpose, the Finite-Difference Time-
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Domain Method (FDTD) was used. TMs method is used to solve Maxwell’s equations
[15][16]. FDTD has proven to be an effective technique of calculating the interaction
of electromagnetic waves with bodies of different material, and complex, hence realistic
geometries [17].
FDTD is considered to be one of the most popular and robust means to
implement engineering electromagnetics models [18].
Beginning with the study of
microwave absorption within the human eye [19], FDTD has been an increasingly
popular tool to study electromagnetic wave interactions with biological tissues for both
impulsive and time-harmonic cases. Such bio-electromagnetics studies have generally
been in the context of partial-body of whole-body simulations of humans exposed to
field sources ranging from power lines to RF and microwave antennas to sources of
visible light.
In this work, a numerical study using the Finite Difference Time Domain
method has been carried out to analyze the use of incident EM fields at microwave
frequencies for the specific case of lung cancer detection. For this purpose, a
computational model of a human torso slice is generated. Several dimensions and
positions for a tumor in the lung area are simulated in order to fmd the optimum
frequency to detect the presence of this malignant tissue. Based on the results presented
in this thesis, 3 GHz appears to be the optimum frequency, in the range 1 - 3 0 GHz, for
detection of lung cancer, according to value of the signal and location of the signal with
respect to the tumor. A time series analysis is also suggested to determine the optimum
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
time, during the incidence of a plane wave on the human torso slice, to detect the
cancerous tissue.
The organization of this thesis is as follows: chapter 2 offers some information
about the lung cancer, such as tumor location, growing process, types, etc. In chapter 3,
a description o f the numerical technique used, namely, the FDTD method is given.
Chapter 4 presents the validation of the accuracy of the program used to perform the
simulation, by comparing the results obtained using our FDTD code with both
analytical and published results.
Chapter 5 explains the approach to generate the
computational model. The simulations, obtained results, and analyses of the results are
presented in chapter 6. Finally, the contribution of this work and suggestions for future
work related to this thesis are discussed in chapter 7.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 2
LUNG CANCER INFORMATION
Lung cancer is known as one of the most difficult cancer to cure, and the
number of deaths that it causes is usually increasing, as the disease is detected at its
evaluate stage. In this chapter, a brief description of the lung cancer is given, including
growing process, types, and location.
Additionally, a brief description of screening
methods is also presented.
2.1
Function of lungs.
The lungs are part of the respiratory system and fill most of the chest cavity.
The organs are irregularly cone-shaped and are separated from each other by the heart,
windpipe and esophagus.
The right lung has three sections and is a little larger than the left lung, which
has two sections. Each normal lung is about the size of a hand and together they weigh
about two pounds. The lungs are filled with air and are covered with a membrane
called the pleura [1].
The function of the lungs is the continuous exchange of gases between the body
and the atmosphere. They alternately cause to exhale carbon dioxide, a waste product
of the cells of the body, and inhale oxygen, a necessity for cellular function.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The
contraction of the diaphragm (the main respiratory muscle) draws air through the nose
and mouth, down the trachea (windpipe) into the lungs where gas exchange occurs. Gas
is brought to one side of the lung by the airways and blood is brought to the other side
by capillary blood vessels. When inhaling, the upper airways immediately begin to
filter the air before it travels down into the lungs. When exhaling, carbon dioxide is
expelled.
f|
f
Figure 2.1 Respiratory system [20].
The lungs, just like all of our organs, are composed of cells. There are several
different kinds of cells found in the lining of the lungs. Some of these cells produce
mucus inside the bronchi and bronchioles (the air ducts that transport air to the lungs)
and trap foreign particles that we breathe into our lungs.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2.2
Description of lung cancer.
Like all cells of the body, the cells in the lungs-divide and reproduce at a
controlled rate to repair worn-out or injured tissues and allow for normal growth. Lung
cancer develops when cells inside the lungs multiply at an uncontrollable rate. These
abnormal tissue masses are called tumors. Tumors are either non-cancerous (benign) or
cancerous (malignant). Figure 2.2 shows: a) the normal distribution of cells in a tissue,
and b) when the cells are forming a tumor.
a)
b)
Figure 2.2 a) Normal cells; and b) ceils forming a tumor [20].
Benign tumors may grow causing discomfort and bleeding, but they do not
spread to other parts of the body and generally are not life-threatening.
Malignant
tumors spread into other areas of the body and destroy normal tissue. As cancer cells
break away from the primary (original) tumor in the lung and metastasize (spread)
through the blood and lymph systems, they form metastatic (secondary) tumors. Lung
cancer can spread to any organ in the body.
Most often secondary tumors will develop in the brain, bone, liver and bone
marrow. Secondary tumors are referred to as metastatic lung cancer rather than brain or
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
iiver cancer to indicate that they are part of a single disease originating in the iungs and
are not new cancers starting in these organs.
2.3
Types of lung cancer.
Lung cancers are divided into two major groups [1], which make up more than
90% of all limg cancer cases; namely, small cell lung cancer (SCLC) and non-small cel!
lung cancer (NSCLC).
Squamous cell, adenocarcinoma and large cell carcinoma are generally grouped
together as "non-small cel! lung cancer" (NSCLC). The way non-small cell lung cancer
spreads and its treatment is very different from that of small cell lung cancer.
Sometimes a lung tumor contains more than one type of cancerous cell.
Because of the different types of lung cancers and different ways they may respond to
treatment, it is very important to identify the cell type so that the best treatment can be
given.
Small cell lung cancer (SCLC) makes up about 20% to 25% of all lung cancer
cases. The identified cells include oat cell, intermediate cell (cells which have many
sides and are spindle shaped) and combined cell (small cell combined with other types
of cancer). Oat cell and intermediate cell types are often found together in the same
tumor.
Small cell lung cancer is most often found in the bronchial submucosa (a layer
of tissue beneath the epithelium, the mucousmembrane) and is more likely to have
metastasized (spread) at the time of diagnosis than other types of lung cancer.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Small cell lung cancer is separated from other cell types due to its rapid growth
rate. Symptoms are of brief duration prior to diagnosis. Eighty percent of small cell
lung cancer is located centrally and 20% is in the periphery of the lung. Histologic
verification (microscopic study of the tissue structure) by a pathologist of small cell
lung cancer is mandatory because treatment is significantly different from nonsmall cell
cancer. The characteristics of small ceil lung cancer are summarized in table 2.1.
Table 2.1 Characteristics of Small Cel! Lung Cancer.
Percentage of occurrence:
20% of all lung cancers
Area located in lung:
Hilar
Progression rate:
Very fast
Inclined to spread:
Very early, (metastases frequently present when
diagnosed)
Epidermoid or squamous cell carcinoma is the most common form of lung
cancer worldwide, accounting for over 30% of lung cancers. Squamous cell carcinoma
usually starts in the large bronchi and very often stays in the chest, without spreading,
for longer periods of time than other lung cancers. The characteristics of squamous cell
lung cancer are summarized in table 2.2.
Table 2.2 Characteristics of Squamous Cell Lung Cancer.
Percentage of occurrence:
30% of all lung cancers
Area located in lung;
Hilar
Progression rate:
Slow
Inclined to spread;
Late, primarily to hilar nodes
10
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
In the United States adenocarcinoma is increasing and surpassing squamous cell
in the number of reported cases. Adenocarcinomas account for 35% of all lung cancers
and have cube or column-shaped cells. These tumors are often found along the outer
edges of the lungs and under the lining of the bronchi. New studies are under way to
examine reasons for the increase in the rate of adenocarcinoma.
In table 2.3 the
characteristics of adenocarcinoma cell lung cancer are summarized.
Table 2.3 Characteristics of Adenocarcinoma.
Percentage of occurrence:
35% of all lung cancers
Area located in lung:
Outer edges of the lung.
Progression rate:
Average
Inclined to spread:
Average.
Large cell carcinoma is a term used for lung cancers that do not fit into other
categories of epidermoid, adenocarcinoma or small cell. Large cell carcinoma makes
up about 15% of lung cancer cases. This type of cancer is most often found in the
smaller bronchi. In table 2.4 the characteristics of large cell lung carcinoma cancer are
summarized.
Less prevalent lung cancers are carcinoid tumor, accounting for 20% of lung
cancers, and bronchioalveolar carcinoma, a subtype of adenocarcinoma.
Carcinoids
arise from glands near the bronchi, while bronchioalveolar carcinoma develops around
scars on the outer edges o f the lungs. Figure 2.3 depicts a CT scan showing a tumor.
11
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table 2.4 Characteristics of Large Ceil Carcinoma.
Percentage of occurrence:
15% of ail lung cancers
Area located in lung:
Central of outer edges of the lung.
Progression rate:
Very fast
Inclined to spread:
Early
Tumor
Figure 2.3 CT scan showing a tumor [21].
2.4
Screening of lung cancer.
Screening is defined as looking for a disease before symptoms appear. Two
tests are usually employed to screen lung for cancer. One test is the routine chest x-ray.
The other main screening test is sputum cytology [1].
Early detection of lung cancer is critical to improving chances of survival.
Physicians use a number of different tests to detect and diagnose lung cancer, including
sophisticated imaging scans that provide more accurate and sensitive results than
conventional X-rays.
The information from these tests enables the physician to
determine the type and stage of the cancer and the best way to treat it [22].
The current tests to perform the lung cancer screening include;
12
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Physical ExamiBation.
This examination is important for detecting any signs of cancer
such as swollen lymph nodes in the neck or collarbone area and also for evaluating
overall state of health.
Chest examination. - Examining the chest and listening to the lungs with a stethoscope
provides information about abnormal breathing sounds or patterns.
Chest X-ray. - X-rays are "flat" pictures of the lungs, which help to identify abnormal
growths [23]. A chest x-ray may not show a tumor, possibly because it is too small or it
may be hidden behind a rib or the breastbone. But, the x-ray may show other clues that
indicate a problem related to lung cancer.
For example, an x-ray may show an
accumulation of fluid between the lung and the chest wall, called a pleural effusion. An
x-ray may show enlarged lymph nodes or pneumonia. New x-ray technologies may
help to make chest x-rays more useful for diagnosing lung cancer. One of the new
technologies allows x-rays to bypass the bones in the chest so that only the lung tissue is
seen, making it easier to read the x-ray.
Computed tomography (CT) scan. - This test is a sophisticated instrument that uses a
computer to create a two-dimensional scan from a series of X-ray images [23]; the
newest version of the CT is called a helical (or spiral) scan. CT scans reveal much more
detail than x-rays and the new helical scans are even more sensitive than regular CT
scans. An improved CT technology is a “spiral” or “helical” CT scan [24].
Positron Emission Tomography (PET) scan. - This test is a scan that traces the way the
body cells act on sugar. PET scans can find cancerous tumors because of their ability to
take up radioactive sugar. PET scanning is a relatively new technology. It is different
13
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
from CT and MRI scanning because PET scan discriminates between cells that are
rapidly dividing, such as tumor ceils, and normal ceils. Currently, research is conducted
to obtain the most information possible from the image [25].
Magnetic Resonance Imaging (MRI). - This technique is similar to a CT scan except it
uses a magnetic field in place of X-rays to create an image [23]. The MR experiment,
in its most basic form, can be analyzed in terms of energy transfer. The patient or
sample is exposed to energy at the correct frequency that will be absorbed. A short time
later, this energy is reemitted at which time it can be detected and processed [26].
Sputum cytology. - In this test a sample of sputum or mucus is obtained after deep
coughing. The sample is then reviewed under a microscope with special stains in order
to determine if cancerous cells are present. While abnormal cells in the sputum did
reveal cancers that were too small to be seen on x-rays, the cells could not reveal the
exact location o f the tumor. Further examination to locate the tumor site is, therefore,
necessary and is often expensive.
Bronchoscopy. - Viewing of the lungs through a hollow, flexible tube (bronchoscope)
that is passed through the nose and throat into the main airway of the lungs.
If
abnormal areas or tumors are seen, biopsies can be obtained through the bronchoscope.
Biopsy. - Removal o f a lung tissue sample for examination under a microscope.
Biopsies are obtained in different ways depending on the location of the tumor: a)
through a bronchoscopy, b) by inserting a needle through the chest into the lung, c) by
removal and examination of an enlarged lymph node in the neck, and d) by a small
surgery on the lung.
14
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
In conclusion, the newer imaging scans (CT, PET and MRI) are very sensitive
and can reveal cancerous growths not seen by conventional chest X-rays. Clinical trials
are underway to determine the effectiveness of screening to permit the early detection
of lung cancer based on these new advances.
2.5
Lung cancer facts.
The following is a list of a few facts related to the lung cancer [22]. Needless to
say, these are reasons to look for new methods to detect the lung cancer in its early
stage.
•
Lung cancer will kill approximately 68,800 w^omen in the U.S. this year more than breast and ovarian cancer combined.
Lung cancer deaths
surpassed breast cancer deaths in 1987.
•
In 1950, lung cancer accounted for only 3% of all cancer deaths among
women; however, by 2000, it accounted for an estimated 25% of cancer
deaths—an estimated 1 of every 4 deaths due to cancer and nearly 1 of
every 8 newly diagnosed cancers among women.
•
Research shows that women are approximately 1.5 times more likely to
develop lung cancer than men.
•
Survival was higher among women with localized disease (52.5%), but
only 16% of cases among women were diagnosed at this early stage.
•
Survival rates declined with age at diagnosis and advanced stages but were
higher among women than among men at all ages and stages and for all
cell types.
15
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
® Recent advances in treatment have increased survival rates and improved
the quaiit>' of life of patients significantly.
•
Thanks to new treatments, the cure rate for lung cancer has doubled over
the last 30 years.
»
Early detection is critical. When caught early, most lung cancers can be
cured.
16
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 3
FDTD DESCRIPTION
It is indispensable to solve Maxwell’s equations in order to study the
electromagnetic fields in different media. One manner to get the solutions is by using
numerical methods. Numerical methods did not gain much importance until recently
because of the lengthy computational time and enormous amount of memory required to
get solutions. But now, with the current powerful computational resources, numerical
methods are becoming the primary methods to observe electromagnetic wave
propagation, scattering and penetration.
The Finite Difference Time Domain method (FDTD) technique has been shown
to be ideally appropriate to compute the electromagnetic response from an object of
arbitrary shape and dielectric value. In this chapter, a description of FDTD is presented.
3.1
Introduction
The Finite Difference Time Domain method has many formulations; scattered
field, total field, potential, implicit, etc. The most common formulations in FDTD are
the total field formulation and the scattered field formulation [27].
While integral equation based methods have been employed occasionally in
time-domain applications, partial differential equation based methods have proven to be
17
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
more practical and have thus seen broader use. Here it is discussed in some detail the
Cartesian mesh finite-difference time-domain (FDTD) jnethod, the oldest and most
widely used of these techniques. This method was first described by K. Yee in 1966
[15]. As computing hardware became faster and more widely available, the method
found numerous applications in electromagnetic radiation, scattering, and coupling.
Stimulated by these successes, various researchers developed enhancements to the
original technique that greatly extended its range of applicability.
X
Figure 3.1 Positions of field components. The E-components are in the middle of edges
and the H-components are in the center of the faces.
Several keys attributes combine to make the FDTD method a useful and
powerful tool. First is the simplicity of the method; Maxwell’s equations in differential
18
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
form are discretized in space and time in a straightforward maimer. Second, whereas a
majority of the numerical electromagnetic modeling techniques are in the frequency
domain, the FDTD method is a time domain solution. Hence, since the method tracks
the time-varying fields throughout a volume of space, FDTD results lend themselves
well to scientific visualization methods. These, in consequence, provide the user with
excellent physical insights on the behavior of electromagnetic fields. Therefore, the
FDTD simulation gives us the opportunity to observe the wave as it propagates through
a medium, interacts with another, penetrates into it and also scatters of it.
This
observation can be done at any point in space and at any time.
Finally, the geometric flexibility of the method permits the solution of a wide
variety of radiation, scattering, and coupling problems.
Yee showed [15] the positions of the various field components in order to meet
the requirements for FDTD. Figure 3.1 shows these positions. In Table 3.1 are listed
the different fields locations in the Yee cell.
Table 3.1 Location of the different fields in the Yee cell
Element
Location
x = (/ + l/2)A x, j = 7 Ay, z = k h z
E % ],k )
x = i Ax, y = { j + l/2)A_y , z = kAz
E ’A h h k )
x-ihx, y -
7
Ay, z = (k + l/2)A z
x = i A x , y = (7 + 1/2)Ay, z = (k + l/2)A z
x = (i+ l/2)A x, y = j A y , z = (k + l/2)Az
x = (z + l/2)A x, y = (7 + 1 /2 )A y , z = kAz
19
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.2
One-dimensional free space formulation
To give a description of the FDTD method, initially a ID-formulation will be
presented. Although the basic idea is the same for different formulations, the approach
used in [28] will be explained here.
The time-dependent Maxwell’s curl equations in free space are
(3.1)
dt
Gn
(3.2)
dt
Ao
E and H are vectors in three dimensions, so in general, equations (3.1) and (3.2)
represent three equations each.
dE^ _ 1
dH,
dH^
dt
dy
dz
dE y
dH,_
dH,
dt
dz
dx
_1_ dH^
dt
dx
^0
dH,
dE,
dz
dy
d ^
J_ (dE,
dE,
dt
Mo ^ dx
dz
J _ f dE
dE y
dt
Mo
Mo
I
^
(3.4)
(3.5)
dy
J_ ( be.
dt
(3.3)
dx
20
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.6)
(3.7)
(3.8)
The explanation will start with a simple one-dimensional case using only E
and Hy , so eqs. (3.3) and (3.7) become
1 dH y
dE^
dt
(3.9)
Sfs dz
=
dt
(3.10)
Pq dz
These equations represent a plane wave with the electric field oriented in the x
direction, the magnetic field oriented in the y direction, and traveling in the z
direction.
As next step, the free space field equations are differenced. In essence finite
differencing replaces derivatives with differences (Appendix A):
Sf _
fixJi)
f{x ,t^)-f{x,t,)
At
dt
At
^ - lim / ( ^ 2»0 - / ( ^ i >0 ,, /(X 2, t ) - / ( x i , t )
dx
Ax
Ax
Taking the central difference approximations for both the temporal and spatial
derivatives gives
E f ‘\ k ) -
1
At
H ;{k + H2) - h ;{k - l l 2 )
Sg
At
Az
Pq
Az
where in the above approximation At and Ax are finite rather than infinitesimal.
21
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.13)
(3.14)
It is used the central difference formula, because of the greater accuracy, since
error decreases as the second power of At and Ax .
In these two equations, time is specified by the superscripts; i.e., “n” actually
means a time i = A t - n . The term “ « + l ” means one time step later. The terms in
parentheses represent distance, i.e., "k" actually means the distance z - A z - k .
The
formulation of eqs. (1.3a) and (1.3b) assumes that E and i f fields are interleaved in
both space and time. Field H uses the arguments A:+ 1/2 and k - H 2 to indicate that
the H fields values are assumed to be located between the E field values. Similarly,
the « + 1/2 and « - 1/2 superscript indicates that it occurs one half time step after or
before n , respectively.
Equations (3.13) and (3.14) can be rearranged in an iterative algorithm:
( k ) ----- — H ’' ( k + l / 2 ) - H " ( k - \ / 2 )
s^-Az
E f ’^ (k) =
— b
//fl •Az
H f {k + l / 2 ) = H l i k + l l 2 )
X
✓
/
/
#
g n -1
f (k + l ) -
(3.15)
{k)\
s
\
tth-1/2
^
I
Trn+1:
E” V ^
\
/
\
\
✓
/
Figure 3.2 Leapfrog time integration.
22
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.16)
Notice that the calculations are interleaved in both space and time. In eq. (3.15),
for example, the new value of
most recent values of
is calculated from the previous value of
and the
. This is the fundamental paradigm of the finite-difference
time domain (FDTD) method [15].
This spatial and temporal interleaving is because of the central differencing.
What results is often referred as "leapfrog" scheme [29] (Fig. 3.2).
Equations (3.15) and (3.16) are very similar, but because
several orders o f magnitude,
and
and
differ by
will differ by several orders of magnitude.
This is circumvented by making the following change of variables:
(3.17)
E =J ^ E
1/^0
This is a system called Gaussian units. The reason for using it here is simplicity
in the formulations. The E field and the H fields have the same order of magnitude.
This has an advantage in formulating the perfectly matched layer (PML), which is a
crucial part of FDTD simulation.
Substituting this into eqs. (3.15) and (3.16) gives
E f ^ {k) =
{k) —
----- — H I
(3.18)
+
1^0-Mo ^
i f f ’ ( i +1 / 2) = i f ; (A:+ 1/2) —
I
—
(^ +1) -
(i)]
23
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.19)
3.3
Determining the cell size
The choice of the cell size is critical in applying FDTD.
It must be small
enough to permit accurate results at the highest frequency of interest, but large enough
to keep resource computational requirements manageable. Cell size is directly affected
by the materials present. The greater the permittivity or conductivity, the shorter the
wavelength at a given frequency and the smaller the cell size required.
A constraint used frequently is ‘TO cells per wavelength”. This means that the
side of each cell should be 1 /1 0 or less at the highest frequency (shortest wavelength)
of interest.
Once the cell size has been determined, the number of cells needed to model the
object and a reasonable amount of free space between the object and the outer boundary
in each dimension is found, and from this the total size of the FDTD space is
determined.
3.4
Time step size for stability
Once the cel! size is defined, the maximum time step At is determined using the
Courant stability condition (eq. 3.20). A larger time step results in instability.
To
understand the basis for the Courant condition, consider a plane wave propagating
through an FDTD grid (Fig. 3.3). In one time step any point on this wave must not pass
through more than one cell, because during one time step FDTD can propagate the
wave only from one cell to its nearest neighbors.
24
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
\
\
Ax
Ax
Plane wave
Figure 3.3 Plane-wave propagation along the diagonal of a two-dimensional
square-cell ( Ax x Ax) mesh is equivalent to one-dimensional propagation
along a mesh with cell size h x ! 4 2 , thus accounting for the
two-dimensional stability limit of c A t < A x / 4 2 [29].
To determine the time step constraint, a plane wave direction is picked so that
the plane wave propagates most rapidly between field point locations. This direction
will be perpendicular to the lattice planes of the FDTD grid. For a grid of dimension d
(where d =1, 2, or 3), with all cell sides equal to A u , it is found that with v the
25
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
maximum velocity of propagation in any medium, usually the speed of light in free
space,
(3.20)
vht < ~
■\jd
for stability.
Depending on the analysis, the Courant condition will have the following
values:
M <—
Cq
for 1-D
(3.21)
At <
for 2-D
(3.22)
for 3-D
(3.23)
-h e.
At <
h e ,
assuming Ax = Ay - A z , which is commonly used.
However, for convenience the
following determination of At is used
=
(3.24)
2 • Cq
After determining the time step, a number of time steps sufficiently large is
estimated to allow a full characterization of the interaction of object and fields.
3.5
Absorbing boundary condition in one dimension
Absorbing boundary conditions are necessary to keep outgoing £ and H fields
from being reflected back into the problem space. Normally, in calculating the E field,
the surrounding H values are needed; this is a fundamental assumption of the FDTD
26
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
method.
At the edge of the problem space the value to one side will not exist.
However, it can be used the fact that there are no sources outside the problem space.
These two facts are used to estimate the value at the end by using the value next to it.
The analysis starts by looking for a boundary condition at the end where k = 0 .
If a wave is going toward a boundary in free space, it is traveling at
, the speed of
light. So in one time step, eq. (3.24), of the FDTD algorithm, it travels
c„-A/ = c , . ^ = ^
2 -Co
2
(3.25)
This equation basically explains that it takes two time steps for a wave front to
cross one cell. Hence, an easy approach to obtain a boundary condition could be
£^"(o) = a ;-2( i )
(3.26)
Although this approach is only appropriate for one-dimensional FDTD
simulation, the same idea is used in 2D- and 3D-FDTD, when an incident plane wave is
generated.
3»6
Propagation in a dielectric medium
In order to be able to simulate a medium with a dielectric constant other than
one, which corresponds to free space, the relative dielectric constant
has to be added
to Maxwell’s equations:
dE
1
dt
■^0
VxH
(3.27)
~ =~— V xE
dt
/i„
(3.28)
27
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Staying with the one-dimensional example and make the change of variables in
eqs. (3.9-3.10) the following equations are obtained
1
BEM
dt
(3.29)
dz
1
dH j4)_
d E . it )
(3.30)
dt
and then go to the finite difference approximations
H l ^ \ k + l / 2 ) = H " (k + l / 2 ) — p = i = - —
3.7
(k + l ) - ( k ) ]
(3.32)
Propagation in lossy dielectric medium
In sections 3.2 and 3.6 the analysis for EM propagation in free space or in
simple media that is specified by the relative dielectric constant
was performed.
However, there are many media that also have a loss term specified by the
conductivity. This loss term results in the attenuation of the propagating wave.
To do this analysis, the time-dependent Maxwell’s curl equations will be written
in a more general form, which will allow to simulate propagation in media that have
conductivity:
8F
s ~ =V xH ~ J
dt
=
1
V xE
dt
J is the current density, which can also be written as
28
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.33)
(3.34)
(3.35)
J = a-E
where a is the conductivity. By putting this into eq. (3.33) and dividing through by the
dielectric constant the following equation is obtained
^
=
(3.36)
— — E
S’o^r
dt
This equation is now reverted to the simple one-dimensional equation:
dE,{l)
1
a
dt
s^Sf^
dz
£^Sq
and by making the change of variable gives
1
dE^{i)
dHy {t )
CT
(3.38)
-
(3.39)
V^o/^o
dt
Next, by taking the finite difference approximations for both the temporal and
spatial derivatives, the following equations are obtained:
At
Az
S^JsnjUn
CT E f ' \ k ) + E 1- ^ ' \ k )
2
g^Sq
In eq. (3.40), some constants can be simplified by using eq. (3.24), hence
1
At
1
Therefore,
29
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.40)
M iJ
1-
At a
1+
(3.40a)
le^e^
1/ 2
1+
3.8
Ata
2 s ^£ q
Reformulation using the flux density
The form of Maxwell’s equations given in eqs. (3.1) and (3.2) has been used,
which uses only the E and H fields. However, a more general form is
3D
(3.41)
■= V x H
dt
D{o}) = £f^ -e* {{(>)■ E { co)
3H
1
,=
dt
„
vxE
jI q
(3.42)
(3.43)
where D is the dielectric flux density. Notice that eq. (3.42) is written in the frequency
domain. The reformulation using the flux density will be started by normalizing these
equations,
n r
E = \^ E
1
D =
^0
(3.44)
-D
(3.45)
■Mo
which leads to
1
dD
rVxH
^oMo
30
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.46)
D{a)) = £l{e})-E{a))
(3.47)
^
(3.48)
=
This form of equations will lead to the very simple finite difference equations.
The only change is the use of D instead of E . However, eq. (3.47) has to be taken into
a time domain difference equation for implementation into FDTD. The first task is to
get it from the frequency domain to the time domain. It is assumed that the medium is a
lossy dielectric of the form
< (» ) = £ ,+ - ^
(3.49)
•E{(o) +
(3.50)
jms^
and substitute eq. (3.49) into (3.47):
D{ q}) =
- •■E { a )
jo)£o
Taking the first term into the time domain is simple because is a multiplication,
hi the second term, knowing that 1/ ja> in the frequency domain is integration in the
time domain, so eq. (3.50) becomes
D{t) = e ^ - E { t ) + — l E { t ' y d t ' .
(3.51)
In order to have an expression in the sampled time domain, the integral is
approximated as a summation over the time steps A t:
D"
- E ” + ^ - ^ ^ E ‘.
^0
(=0
31
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.52)
Note that E and D are specified at time t = n-A t. Since it is needed to solve
for
given the value o f D \ and E ” is needed in the calculation of the summation,
the E ” term is separated from the rest of the summation:
(3.53)
D ‘^ = s ^ - E ’'
^0
^0
!=0
Now E" can be calculated from
E'
--------- — .
<7-At
s ^ + -------
(3.54)
^
The current value o f E , E " , can be calculated from the current value of D and
previous values of E . For convenience a new parameter for the summation is defined
/ ■=—
(3. 55)
so eq. (3.54) can be reformulated with the following two equations:
jyn _ jn - \
E" = — ---- i - —
a-A t
(3.56)
^
£, + -------
rn rn-\ O'• At
r= r-^+
E\
’
(3.57)
®’o
Note that the summation is calculated by eq. (3.56), which, at every time step n,
simply adds the value E ” times the constant term to the previous values of the
summation at n-1. The entire FDTD formulation is
dx^=dxj^+Q.5* ( % _ i- % )
* {dXk -
)
32
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.58)
(3.59)
= ix^. + gbxj, * eXj^
h-Vk = %* + 0.5 * {eXf, - ex^^j )
(3.60)
(3.61)
where
gbx^ = a * h t I
(3.63)
Hence, all of the information regarding the media is contained in eqs. (3.62) and
(3.63). For free space, gax = 1 and gbx = 0; for lossy material, gax and gbx are only
values of dx^ and previous values of
3.9
in time domain.
FD ID in two dimensions
The FDTD method can be expressed dimensionally for different cases (ID, 2D
and 3D) according to the analysis and structure of interest. In this study, a 2D analysis
is conducted. Hence, in this section an explanation of the 2D-FDTD is given.
The analysis starts with the normalized Maxwell’s equations used in section 3.2;
*
V x //
(3.64)
D{ o)) = s*{co)-E{a?)
(3.65)
dH
1 „ p
---- = — ------ v x E
(3.66)
St
33
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
yo + (j-^l)Ay
yo+JAy
Xo
+ iAx
Xo
+ (i+l)Ax
Figure 3.4 The FDTD unit cell for the two-dimensional TM case.
When doing the three dimensional simulation, it is needed to deal with six
different fields:
, E^, E^,
, H ^ , and
. In doing two-dimensional simulation,
it is chosen between one of two groups of three vectors each: ( 1) the transverse
magnetic (TM) mode, which is composed of E^,
, and
showed in Fig. 3.4, or
(2) the transverse electric (TE) mode, which is composed of E^, E ^ , and H , , showed
in Fig. 3.5. Equations (3.64-3.66) are now reduced to
1
dt
dEy
dH^
dx
dy
D{ o} ) = £*{6})' E { a )
dt
(3.67)
(3.68)
8y
(3.70)
34
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
yo + (j+l)Ay-
yo+jAy
-^x
Xo + (i+l)Ax
Xo + iAx
Figure 3.5 The FDTD unit cell for the two-dimensional TE case
It is important that there is a systematic interleaving of the fields to be
calculated. Putting eqs. (3.70), (3.72), and (3.73) into the finite differencing scheme
results in the following difference equations:
At
<
g ; ( i+ i/2 .; ) - g ; ( i- i/2 ,y y
y/^o/^o
Ax
1
(3.71)
H ^ jiJ +1 / 2 ) - - 1 / 2 ) '
Ay
f E :^^'% j+ l)-E f'% j))
i 7 f ( i , i + l / 2 ) - 7 f " ( i ,y + i/2 ) ^ ___ 1
-\/^o/^o V.
At
iff
+ l / 2 , i ) - H ; {i + 1/2, j ) _
1
(
^V
At
(3.72)
AF
(3.73)
^
Using the same type of manipulation as in 1-D case, including
Ax
At = 2-Ca
the following equations are obtained
35
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
d^ij = dZij + 0.5 *
- hx^j +
)
(3.74)
ez,. = g az .j * idz,j - iz,j j
(3.75)
K j = K j + Sbz^j * ez,j
(3.76)
^^i,j = H ,y + 0-5 *
)
(3.77)
= hy,j + 0.5 * [ez,^,j - ez,. )
(3.78)
3.10 Total/Scattered field formulation
The simulation of plane waves is often of interest in computational
electromagnetics. Many problems deal with plane waves. Furthermore, after a distance
on the order of tens of wavelengths, the field from most antennas can be approximated
as a plane wave.
In order to simulate a plane wave in 2D-FDTD, the problem space will be
divided up into two regions, the total field region, and the scattered field region (Fig.
3.6). There are two primary reasons for doing this: (1) The propagating plane wave
should not interact with the absorbing boundary conditions; (2) the load on the
absorbing boundary conditions should be minimized. These boundary conditions are
not perfect, i.e.; a certain portion of the impinging wave is reflected back into the
problem space.
By subtracting the incident field, the amount of the radiating field
hitting the boundary is minimized, thereby reducing the amount of error.
36
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
PML
One-dimensional
incident array
Scattered field region
Incident plane
wave is
substracted out
here.
ABCs
jb
I
Total field region
ja
Source
point
Va
i
II
i
Incident plane
wave is generated
here.
ib
la
Figure 3.6 Total field/scattered field of the two-dimensional problem space.
Figure 3.6 shows how this is accomplished.
There is an auxiliary one­
dimensional buffer called the incident array. Because this is one-dimensionai array, it is
easy to generate a plane wave: a source point is chosen and the incident field is just
added at that point. Then a plane wave propagates away in both directions. Since it is a
one-dimensional array, the boundary conditions are perfect.
In the two-dimensional field region every point in the problem space is either in
the total field or it is not; no point lies on the border. Therefore, if a point is in the total
field but it uses points outside to calculate the spatial derivatives when updating its
value, it must be modified. The same is true of point lying just outside that uses points
37
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
inside the total field. This is the reason for the incident array; it contains the needed
values to make these modifications.
Total field
Ja+l
o
Dz
t
Hy
o
Dz
o
Dz
!
Hy
o
Dz
o
Dz
t
Hy
Hx
o
Dz
1
Dz
Hx
t
Hy
Hx
Hx
ja-1
O
Hx
Hx
ja
t
Hy
Hy
o
Dz
t
Hy
Hx
t
Hy
o
Dz
Hx
t
Hy
Hx
Scattered field
ia+1
ia-1
Figure 3.7 Total field/scattered field of the two-dimensional problem space.
There are three places that must be modified:
1. The Dz value at j = j a or j = j b :
A (h j a) = D, (i, j a) +0.5*
{ j a - 1 /2 )
(3.79)
D, (i, j b ) = D, (i, j b ) - 0.5 * D ,
{jb + 1/2)
(3.80)
38
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
2. The
field just outside j = j a and j = j b :
HXhja-M2)^H,(i,ja-\n)+<i.S*E^j„(ja)
(3.81)
H , ( i J b + \ ! 2 ) = H j 4 , j b + H 2 ) - G . 5 ‘' E , j ^ ( j b )
(3.82)
3. The H t field just outside i = ia and i = ib :
/ f , ( i a - i / 2 ,; ) = / / , ( f a - l / 2 ,y ) - 0 . 5 * £ ._ , „ ( i )
(3.83)
O ')
ff (/6 +1 / 2. y) = H (ib +1 / 2, y) + 0.5 •
(3.84)
Figure 3.8 shows the results of a simulation of a pulse plane wave propagating
along the total-field FDTD. This pulse is generated at the boundary of the scattered/total-field space in one end, and subtracted at the boundary on the other end. The
shown sequence was recorded at 30, 60, 90 and 115 iterations.
T*80
S3 60
S3 60
cm
cm
T = 1 1 5 “ ” T..........
0,5 T
^
23
40
cm
m
.
23
BO cm
■
20
40
^
cm
60 60
cm
Figure 3.8 2-D simulation of a plane wave propagating in free space. The pulse is
generated at one end and subtracted out at the other end.
39
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
3.11 Absorbing boundary conditions for 2D case
To use the FDTD method to solve open-region problems, like modeling a
stmcture that is situated in free space, such as a scatterer or radiating antenna, it is
required some mechanism for truncating the mesh at a finite distance from the
computational region of interest. Ideally, this mechanism, namely absorbing boundary
condition, should permit electromagnetic energy to pass out of the problem space
without distorting the fields or reflecting energy back into computational domain.
An absorbing boundary condition may not always be necessary when applying
FDTD. If the FDTD problem space is bounded by a condition that can be implemented
directly into the finite difference equations, an ABC is not necessary.
Another consideration is that of determining the distance between the object and
the outer boundary. The farther from the object the outer boundary is located the better
the absorption of the outward traveling waves. This is due to the fact that these waves
become more like plane waves as they travel farther from the structure that radiates
them. However, the number of cells that can be placed between the object and the outer
boundary is limited by computer memory.
There have been numerous approaches to this problem [27,30]. Among these
approaches, one of the most flexible and efficient ABC’s is the perfectly matched layer
(PML) developed by Berenger [31], which employed a fictitious, directionally
dependent pair of electric and magnetic conductivities for the purpose of absorbing
outgoing waves and minimizing the reflection back into the problem space.
40
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The approach used in the simulations performed in this thesis is a slight
modification to the original Berenger implementation. This modifieation was proposed
by Sullivan [32]. It is different in two ways: first, it utilizes a fictitious conductivity
associated with the electric displacement, instead of the electric field. The purpose for
using displacement is that those conductivities connected with the PML are completely
separate from any electrical conductivity associated with the problem space. Secondly,
it suggests varying the FDTD parameters directly as the most efficient way of
implementing the PML, as opposed to varying the conductivities or the size of the cells,
and then converting it to FDTD parameters.
The basic idea for the PML approach is the fact that if a wave is propagating in
medium A and it impinges upon medium B, the amount of reflection is dictated by the
intrinsic impedances of the two media
r =^ -llg .
(3.85)
nA+VB
which are determined by the dielectric constants e and the permeability pi of the two
media
ri = ^
(3.86)
In many situations it is assumed that // is a constant, so when a propagating
wave goes from
to
, it meets a change in the impedance, therefore a portion of the
wave is reflected. This proportion is given by eq. (3.85). However, if pi changed with
£ so J] remained a constant, F would be zero and no reflection would occur. This
solves only one part of the problem, because the wave will continue propagating in the
41
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
new medium. This new medium must be also lossy so the wave will fade away before
it hits the boundary. This is accomplished by making both s and p. of eq. (3.86)
complex, because the imaginary part represents the part that causes decay.
By moving eqs. (3.67) to (3.70) to the Fourier domain, it is obtained
d
dH^
dH,
dx
dy
(3.87)
(3.88)
X&) = s l ( a ) - E,{a>)
j( o H , = -C(
ja H
DE.
(3.89)
(3.90)
dx
Note that s and p were eliminated from the spatial derivatives in eqs. (3.87),
(3.89), and (3.90) for the normalized units. Instead of putting them back to implement
the PML, fictitious dielectric constants and permeabilities
jQ)D,-Sp^{x)-ep^iy) = c,
, p \ , and p*p^ are added:
dH^
dH,
dx
dy
(3.91)
(3.92)
jo )H , ■p p , (x) ■pp, (y) = -Co
dE.
dy
jo)Hy ■Ppyix)-ply{y) = c^
dx
(3.93)
(3.94)
Here is worth to note that the value Sp is associated with the flux density D ,
not the electric field E ; and, two values have been added each of Sp in eq. (3.91), and
Pp in eqs. (3.93) and (3.94), one for the x direction and one for the y direction; and
42
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
finaily, nothing was added to eq. (3.92). These fictitious values to implement the PML
have nothing to do with the real values of s'l {(o) which specify the medium.
After several steps, the equations to represent D , , H^,
o r ' " (h j ) = gl3(i) •gj3{j) ■
are the following
j)
+ gi2(i) ■gj2{j) ■0.5 ■
(3.95)
[A;(i + i / 2 , j ) - i f ; ( f - i / 2 , / ) - i f ; ( u - + i / 2 ) + i / ; ( i , j - i / 2 ) ]
i T f *(i + 1/2, j ) = fi3{i +1 / 2)- i f ; (/ +1 / 2, j ) - fi2{i +1 / 2) •0.5 •curl _ e
(3.96)
+#(i)-/";‘''^(i + l/2 ,i)
where
(i +1 / 2, i ) = r~l'^ 0 + 1/ 2, j ) + curl _ e
curl_e =
+ 1)
(3.97)
(3.98)
and
i f f (i,i + l/2 )= ii3 ( i + l / 2 ) - i / ; ( i , i + l /2 )
+ j j 2 { j + l/2)-0.5-CM r/_e
(3.99)
+yii(0-/r^(i,i+i/2)
where
' (i, 7 +1 / 2) = / f (i, J +1 / 2) + cmH _ e
(3.100)
The values of the M I set of parameters associated with the PML as well as a
complete description of the steps to apply the simplified PML are given in [28] [32].
Figure 3.9 shows the results of a simulation to illustrate qualitatively the
effectiveness of an 8-point PML with the source offset from the center in both the v and y ~ directions. Note that the outgoing contours remain concentric. Only when the
wave gets within the PML, the distortion starts to take place. The efficacy of the PML
43
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
is evident in the bottom figure because the contours would not be concentric circles if
the outgoing wave were partially reflected.
In chapter 4, the effectiveness of the absorbing boundary condition is shown
quantitatively when comparing FDTD results with analytical solutions.
SO
50
40
UJ
30
20
-0.5 i :
10
21
cm
40
m
cm
T *1®
E 3 0 j..
Figure 3.9 Results o f a simulation to show qualitatively the effectiveness of the
absorbing boundary condition.
3.12 Numerical dispersion
The numerical algorithms for Maxwell’s curl equations defined by the finitedifference systems cause dispersion of the simulated wave modes in the computational
lattice. That is, the phase velocity of numerical wave modes in the FDTD grid can
differ from the vacuum speed of light c , in fact varying with the modal wavelength, the
44
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
direction o f propagation in the grid, and the grid discretization. Numerical dispersion is
a factor in FDTD modeling that must be accounted for to understand its operation and
its accuracy limits, especially for electrically large stractures.
Taflove derived the dispersion relation for the second-order Yee fmitedifference time-domain grid [30].
To facilitate this examination of the dispersion relation, initially the analysis is
restricted to one dimension (which corresponds to grid-aligned propagation in higher
dimensions), but similar arguments hold for nongrid-aligned propagation in two and
three dimensions.
In one dimension, the angular frequency and the wavenumber are related as
sm
^mAt'^
^cAt^
2
\Az)
\
^
sin^
f
(3.101)
V
which is the dispersion relation.
By defining the stability the factor as
s=
^cAt^
(3.102)
V AZ y
the dispersion relation becomes
sm
^coAt^
= 5-sm
^ kAz'^
(3.103)
which reduces in the case that A z « X (being X the wavelength) to the ordinary
relation
a) = c - k
that is, the nondispersive characteristic.
45
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.104)
s = 0 .6
■ s = 0 .5
- s = 0.1
D.9
>I
0 .7
0.6
> 0 .5
0.4
0.3
0.2
0,1
ID
IS
2D
25
Vti
Figure 3.10 Phase speed (relative to c) versus number of mesh points per
wavelength for three different stability factors.
Using the dispersion relation, eq. 3.101, both the phase and group velocities can
be obtained [33]. The phase speed is given by
_(o _ 2c
sm
k skAz
/
kAz
5-sin
\
V
(3.105)
2. j j
while the group speed is given by
c
^ kAz'^
—cos
s
da
v„ =■
dk
1
.
-v -sm
2
(3.106)
kAz
J
Figs. 3.10 and 3.11 show the phase and group speed versus the number of points
of the space mesh per wavelength for three different stability factors. It can be observed
that even though those velocities differ from c , the condition of stability is verified, as
long as the spatial resolution is lowered.
46
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
To quantitatively assess the dependence of numerical dispersion upon FDTD
grid discretization in the two-dimensional case, it is taken as an example the twodimensional TM mode, assuming for simplicity square unit ceils ( Ax = Ay = A ) and
wave
propagation
at
an
angle
with
a
respect
to
the
positive
x-axis
{ k = k cos a ; k = k s in a ). Then, the numerical dispersion relation simplifies to [30]
\
T ■ \
• 2^ A - k c o s a \
" 2fA k s in a
= sm
+ sm
2
y
K
2
J
V
/
A V .
- sm
cAt)
(3.107)
Equation (3.107) can be conveniently solved for the numerical wavevector k at
any wave propagation angle a by applying the following Newton’s method iterative
procedure:
~
_ p
sin ^
) + sin ^ (m ,. ) - C
A s k i^ A k i j+ B s in (2 M j. j
s = Q.5
s = 0,1
0.3
O.B
0,7
O.B
0.4
0.3
0.2
Figure 3.11 Group speed (relative to c) versus number of mesh points per
wavelength for three different stability factors.
47
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.108)
where
is the improved estimate of k , k,- is the previous estimate of k , and
A , B , and C are coefficients determined by the wave propagation angle and FDTD
grid discretization:
.
A -cosa
2
A = -----------,
„
B =
A -sin a
,
2
^
C = — sm ----{cAtj
{ 2 J
n ina\
{3.109}
A very good value for k^ to start the iterative procedure is simply 2 n , the wave
number of the correspondmg mode in free space. For this case, it is easily shown that
the numerical phase velocity
is given by
V
=
^
Ort
(3.110)
^ fin a l
where k^^^j is the final result of the Newton’s method iterations.
Figure 3.12 shows the results obtained using this procedure that illustrate the
variation of the numerical phase velocity with propagation angle in a two-dimensional
FDTD grid. In Fig. 3.12, three different grid resolutions R of the propagating wave are
examined: coarse (i? = 5 cells/wavelength), coarse-normal {R - 10 cells/wavelength),
and fine-normal {R = 20 cells/wavelength). For each resolution, the time-step relation
cAt = A /2
is maintained.
This relation is commonly used in two- and three-
dimensional FDTD codes to satisfy the numerical stability criterion with ample safety
margin.
48
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0.99
0.98
>
i
0.97
a.
i
0.9B
I
-z.
0.95
0.94
Wave angle,
a (degrees)
Figure 3.12 Variation of the numerical phase velocity with wave propagation angle
in a two-dimensional FDTD grid for three grid resolutions.
From the figure, it can be seen that the numerical phase velocity is always less
than c, the free-space speed of light, reaching a maximum at 45“ (oblique incidence) and
a minimum at 0° and 90° (incidence along either Cartesian grid axis) for all grid
resolutions. This represents a numerical phase velocity anisotropy that is
in h e r e n t
in the
Yee algorithm.
It is also known [30] that a pulse distortion can be bounded by obtaining the
Fourier spatial-frequency spectrum of the desired pulse, and selecting a grid cell size so
that the principal spectral components are resolved with at least 10 cells per wavelength.
This would limit the spread of numerical phase velocities.of the principal spectral
components to less than 1%, regardless of the wave propagation angle in the grid.
49
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 4
COM PAM SON OF FDTB RESULTS W ITH THEORETICAL
AND PUBLISHED RESULTS
In order to test the FDTD program developed for this study, several tests were
carried on. The results were compared with analytical results as well as some published
results.
These tests were; a) Calculation of the reflection coefficient of a dieletric
object, b) Change of value of wavelength of a sinusoidal plane wave when passing from
free space to a dielectric medium, c) Transmitted wave through a dielectric cylinder,
and d) Propagation in biological tissues, namely a human head.
This offers an
opportunity to check the accuracy of the FDTD simulation program used here.
4.1 Calciilation of reflection coefficient
The ratio o f the amplitudes of the reflected and incident electric fields is called
the reflection coefficient and designated by F ,
(4.1)
K
h + 'li
where
/jj = Intrinsic impedance of region 1.
50
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
T]2 - Intrinsic impedance of region 2.
For a perfect dielectric material the intrinsic impedance is
(4.2)
Scattered field region
PML
=4
fx = 0
ii
i
Total field region
Incident plane wave
Free
space
i
Figure 4.1 Configuration o f FDTD space to verify the reflection coefficient for an EM
plane wave traveling from free space to a dielectric medium with Sr = 4 and c = 0.
If region 1 is free space, then
(4.3)
If region 2 has the same permeability but a relative permittivity s^. - 4 , then
Po _ 1 y“o _ 1
„
_
^2
~ if
V
~
O II
2 y £■(,
“
o ^0
2
Therefore the value of the reflection coefficient is given by
51
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.4)
I
r =
(4.5)
'3
Regions 1 and 2 are as shown in Fig. 4.1. The simulation is performed using a
pulse incident plane wave with magnitude 1.
Free
E
Ew
LU
0.5
0 = 0
0
-0.5
1
E
E
T = 75
Q.5
0
-0.5
10
15
20
25
30
35
40
1
E
E
45
50
T = 100
0.5
0
-0.5
10
15
20
25
30
cm
35
40
45
50
Figure 4.2 EM pulse traveling from free space to a dielectric medium with
= 4. As expected, the amplitude of the transmitted wave is 0.66,
while the reflected wave is 0.33.
After 50 iterations, the pulse wave is still propagating only in the free space
towards the dielectric material. Once 75 iterations have elapsed, the pulse has already
hit the boundary between the dielectric material and free space. Finally, after 100
iterations the pulse has been partially reflected and transmitted into the dielectric
52
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
material. As expected, the reflected wave is 0.333 of the amplitude of the incident field
with a phase shift of ^ , while the transmitted wave into the region 2 is 0.666.
4.2 Variation of wavelength.
For this test, the same configuration as in case 1 is used.
However, the
simulation is performed using a sinusoidal incident plane wave with magnitude 1 and
frequency of 600 M H z .
The corresponding wavelength in free space is given by
C(,
3x10^m / s
„
A = — = ---------- r - r - = 0.50 m
®
f
(4.6)
6 0 0 x 1 0 « 5 -‘
The wavelength in a lossless medium can be obtained using the following
expression [34]:
i =^
p
=^
=^
(OPflS
=
fPflS
=
fPflRSR
(4.7)
sJMr S r
For the case of study the wavelength is
A=
Q.5m
^
„
(4.8)
= 0.25 m
Vl-4
In figure 4.3 it can be seen that after 190 iterations, the incident plane wave has
not reached the medium. The wave is still propagating in free space with a wavelength
of 0.50 m .
Once 450 iterations have elapsed, the wave has already passed through the
boundary between free space and the medium with s ^ = A.
wavelength inside this medium is 0.25 m .
As expected, the
As in case 1, the amplitude of the
53
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
transmitted wave is 0.66. ^Tiat is seen in the FDTD space corresponding to free space
is the combmation of both the incident and reflected waves.
1.5
T = 190
1
E
S
\
0.5
e, = 4
0=0
Free sp a c e
^
0
-0.5
-1
-1.5
= 0.5 m
J __________ I
0
20
40
BD
80
100
120
cm
140
1B0
180
200
1BQ
180
20D
T=450
1 = D.25 m
0
20
40
BD
BO
100
120
cm
140
Figure 4.3 EM sinusoidal wave traveling from free space to a dielectric medium with
Sj^ = 4 ,
= 1, and a = O S / m. The wavelength inside the medium changes
accordingly to the permittivity.
4.3 Transmitted wave through a dielectric cylinder.
In this third test for the FDTD program, a plane wave interacting with an object
is simulated, shown in Fig. 4.4. This test is reproduced from [28]. Without loss of
generality, the object structure is specified as dielectric cylinder. The reason for using a
dielectric cylinder is that there exists an analytical solution to this problem. The field
54
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
resulted from a plane wave at a single frequency interacting with a dielectric cylinder
can be calculated through a Bessel function expansion [35] [36].
Scattered field region
PML
Total field region
Central
axis
ii
Incident plane wave
. t . ±
. ±
±
^ ±
Figure 4.4 Configuration of the simulation of a plane wave impinging
on the dielectric cylinder.
The signal to be compared is the amplitude of the electric field
along the
central axis, as indicated in Fig. 4.4.
For a T A f uniform plane wave traveling in the +x direction in free space that is
incident normally on a dielectric circular cylinder of radius a (Fig. 4.5), the incident,
scattered, and transmitted (into the cylinder) electric fields can be written as [35]:
(4.9)
/3=-00
4-co
E: = a,E,
(4.10)
55
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
=a.EoTrc,JAKp)e'jn<?
a.
(4.11)
i / f ^ ( ^ o a ) ( A r ^ a ) / s ^k ^aJ „ (k^ a) - Hf'^' { k ^ a ) l { K ^ ) \
(4.12)
(4.13)
PlSfl©
Incident
wave
-►X
— ®
Figure 4.5 Uniform plane wave incident on a dielectric circular cylinder [36].
The resulting transmitted electric field
from the simulation is going to be
compared to the analytical solution using a Bessel expansion.
The phase constant k
for a perfect dielectric of permittivity' s
and
permeability fi is given by [34]:
■CDJjUS
'
Hence,
is given by
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.14)
(4.15)
0^0
The electrical conductivity of the dielectric scatterer can be accounted by
redefining the permittivity in terms of an effective permittivity as [37]:
(4.16)
l-J-
Thus, the propagation constant for the lossy dielectric can be calculated as
(4.17)
=(o4M dsf
0
20
£
V'
*4:,. ^'4- X
■1 = 2S
«
1 = 80
40
S3**®
30
40
60
S3
40
30
60
cm
0
33
1=75
40
SCh
I
20 mP-tW
1 = 11^
^ * 1?f
20
40
K1
m
cm
40
m
a
cm
Figure 4.6 Simulation of a uniform plane wave incident on a dielectric
circular cylinder.
57
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
For this simulation a material with relative dielectric constant of 30, and
conductivity cr = 0.3 5 / m was selected.
In Fig 4.6 it can be seen the simulation of a plane wave pulse impinging on a
dielectric cylinder. The cylinder has a 20 cm diameter. After 25 iterations the pulse has
begun to propagate. The pulse has already hit the cylinder after 50 steps. Once 75
iterations have elapsed the pulse is partially reflected and transmitted through the
cylinder. Part of the pulse is propagating in free space around the cylinder. After 100
iterations the main part of the pulse has left the FDTD space, but part of it is still
propagating through the cylinder.
O
+
50 MHz
FDTD
Bessel
-10
300 MHz
0.5
-10
700 MHz
-10
-2
0
10
cm
Figure 4.7 Comparison of the FDTD results vs. Bessel function expansion
results along the propagation center axis of a cylinder at three frequencies.
The cylinder is 20 cm in diameter and has a dielectric constant of 30 and a
conductivity of 0.3 S/m.
58
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The results of the simulation for three values of frequencies are compared with
the analytical solution using Bessel fimctions in Fig. 4.7.
The agreement between the theoretical solution and the results of the FDTD
simulation is acceptable.
4.4 Propagation in biological tissues.
Understanding the interaction of EM radiation with humans is essential in a
number of contemporary applications.
There is currently strong interest in the
biological effects and in the medical use of both radio-frequency (RF) and microwave
radiation. Special attention is paid to the absorption of EM energy by the human head.
The use of handheld transceivers for wireless communications, which operate in close
proximity to the head, has raised safety-related questions and questions concerning the
effect of the head on the performance of the mobile phone antenna. Research on the
induced EM field in the human head is aimed at determining the total absorbed power
and the SAR distribution rate.
As a fourth test for the reliability of the FDTD program, the results of a
simulation were compared with the results of the article: “A fmite-element procedure
based on a boundary-value approach for the evaluation of the electromagnetic exposure
in biological phantoms” [38]. The purpose of the simulation is to analyze the electric
field distribution in a human head due to the radiation of an electric field line.
Hence, the first step was to generate the structure to analyze, namely, a
discretized version of a horizontal slice of a human head, located in the vacuum. The
59
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
picture o f the biological phantom (Fig. 4.8) was used to obtain the electrical properties,
namely permittivity and conductivity, of the tissues to analyze.
A y
Tissue properties
O
a
H
B
a
9
Eleclicmagnetic
so u rc ^
H
■
■
■
s = 4 1 .4 - 0 = 0.867 S/m
E = 12.4 - 0 = 0.143 S/m
e, = 42.fi - 0 = 0.782 S/m
s = 5 5 . 2 - 0 = 1 .3 9 4 S/m
^'=68.9 - a = l.fi3fiS/m
e - 55.2 - o = 1.167 S/m
e , = 55.0 - 0 = 0.943 S/m
^ = 61.3 - o = 1.53BS/m
s = 3 2 .5 - a = 0.574 S/m
e - 5 2 .7 - 0 = 0.942 S/m
Figure 4.8 Mode! and characteristics of the tissues in a human head. Figure from
“A finite-element procedure based on a boundary-value approach for the evaluation
of the electromagnetic exposure in biological phantoms” S. Caorsi, E. Bermani and
A. Massa., IEEE Transactions on Microwave Theory and Techniques, Vol. 50, No.
10, October 2002.
After generating the array containing the electrical properties of the scatterer, a
test in this structure was performed using as source an infinite electric line, with
frequency 700 MHz. The simulation of this kind of source is more appropriated to
analyze later the effects of EM interaction of handset antennas and a human being in
personal communications.
60
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
'V
(a)
(b)
Figure 4.9 Electric field distribution |Ez| in the biological structure according to
(a) Finite-EIement Procedure [38] and (b) FDTD.
For purposes of comparison Fig. 4.9 shows patterns of the obtained FDTD
results with the values using Finite-EIement method. A good agreement can be
observed in the whole biological structure. It was not possible to make a quantitative
comparison, due to the lack of more accurate data of the test using the FE method.
In conclusion, the different tests that were performed served as a good reference
to verify the accuracy of the FDTD program used in this work. As they revealed, the
program is reliable.
61
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 5
GENERATION OF THE COMPUTATIONAL MODEL
5.1 Introduction.
The FDTD numerical method is based on a grid [15]. All fields and materia!
compositions appear only at grid points.
Therefore, the scattering object has to be
discretized. This discretization is achieved by first creating the scattering object on the
computer in its real shape, size and material composition and then mapping it onto the
grid. As with most numerical techniques, it is an issue to generate the real object on the
computer.
When working with FDTD, due to the discretization of the object, a
characteristic called staircasing raises. It is therefore important to select an appropriate
cell size to obtain acceptable results, but keeping the computational memory resources
manageable. For illustrative purposes, Fig. 5.1 shows the discretization of a surface. A
not appropriate cell size was used to make the staircasing problem more evident.
For simple objects like cubes, cylinders, etc., object generation is easy since
mathematical formulas are followed. This approach to modeling a solid object, a sphere
for example, is to determine if a particular grid point is within the object being built. If
that is the case, then the array containing the material information is set to the
appropriate values of permittivity and conductivity of the object. A good point to keep
62
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
in mind is that previously to building the object, it is convenient to set the FDTD space
with the electric properties of free space, and then overwrite the computational array,
also called material array in this work, in the determined points according to the object
being built. This w'ay it is ensured that there are not grid points with undetermined
properties.
Discretized object
Continuous object
Mesh
Figure 5.1 Discretization of a surface. For illustrative purposes the ceil size was not
small enough to make staircasing problem more evident.
Contrary to following a mathematical formula, in case of having a complex
object as scatterer, it is a tedious assignment to generate an algorithm, which describes
the object. In some cases, it is possible to superimpose several easy-to-build shapes in
63
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
order to represent a more complex object. However, this is not the case of study of this
thesis, as it will be seen in next section. Hence a different approach has to be used to
generate the object.
5.2 Torso model for sinmlation.
Since the purpose of this thesis is to analyze the use of microwave for detection
of malignant tissue in the lung area, a slice of the human torso had to be generated.
Figure 5.2 shows the image of a slice of a human thorax. The image was
obtained from Stanford University Medical Media & Information Technologies web
page [39]. The objective now is to generate the computational array in our program,
containing the properties of the scatterer, in this case, the thorax, also named torso in
this work.
Figure 5.2 Image of a cross section of thorax from Stanford University Medical Media
and Information Technologies Web Page.
64
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
As an initial stage in this work, the tissues within the torso slice are assumed to
be homogeneous, which means, for instance, that the area representing the lungs has the
same permittivity and conductivity in every grid point.
Furthermore, a simplified version of the torso slice was analyzed. Figure 5.3
shows this simplified model, including the organs within. Similar models have been
used in [40] [41], however in none of them FDTD was used. In fact, they worked with
FEM to analyze current density distribution. Here it is worth to say that there is not
much work reported in the lung cancer detection area using FDTD.
Figure 5.3 Simplified version of a thorax cross section. These organs within the thorax
are: 1) lungs; 2) blood, which fills the atria and the ventricles; 3) bones (ribs and
backbone); 4) heart wall; 5) interstitial fluid; 6) esophagus; 7) intercostal muscles;
and 8) thoracic wall, which includes the skin, the fat, and the muscles.
65
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The approach applied to generate the computational model of the torso was to
analyze the image of the torso. A software application was developed to interpret the
image given in Fig. 5.3. The image was color segmented according to the distinct
tissues represented. Essentially, each pixel in the image represents a point in the FDTD
grid. Since every tissue was related to a different color, then by reading the picture
pixel by pixel, the properties array could be generated.
Original
color
segmented
image.
Rescaling
According to
frequency.
Scanning of
image and
FDTD
simuiatic
generation o f
ID’S file.
Figure 5.4 Flow diagram of generation of the computational torso model.
Figure 5.4 shows the flow diagram of the approach to generate the
computational model. In the next two sections it is described extensively how this
process is performed.
5.3 Computing the torso dimensions in cells.
The purpose of the first software application developed in this chapter was to
compute the dimensions of the torso in cells [Appendix C]. This software application
finds the recommended maximum cell size to be used to build the structure according to
the frequency of analysis. Hence, in order to build the simulated torso, whose physical
dimensions are 38 cm. x 26 cm. (Fig. 5.5), the required number of cells to build this
structure changed according to the frequency of interest.
66
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Given that the higher the frequency, the lower the wavelength, and in
consequence smaller the cell size, it is going to be needed more cells to represent the
same physical m dth and height for a high-frequency analysis, when compared with a
low-frequency analysis.
38 cm.
26 cm,
Figure 5.5 Dimensions of the simulated thorax slice.
Initially, the software application reads the properties, namely permittivity and
conductivity, o f the distinct tissues within the torso, for the range 1-30 GHz.
For
illustrative purpose, Fig. 5.6 shows the plot of those properties for the limgs.
These values were obtained from the Web application from the Institute for
Applied Physics [42].
This web application is aimed to calculate the dielectric
properties of human body tissues in the frequency range from 10 Hz to 100 GHz using
the parametric model and the parameter values developed by Gabriel et al [43 - 45].
67
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
-© - R elati« permittivity
Conductivity fS/m)
20 ^
Frequency (GHz)
Figure 5.6 Lungs dielectric properties for the frequency range 1-30 GHz.
Values taken from the Institute for Applied Physics, Italy; based on
parametric model developed by Gabriel et al[43 - 45].
For purpose of simplicity, as well as because some dielectric properties values
for certain tissues were not available, some tissues were generated by composing them
from other tissues. Obviously, the values of the properties for these composed tissues
were not available directly; hence, they were obtained by making a linear combination
of the components.
The first tissue to be generated using this approach was the thoracic wall, whose
components are skin, fat, and muscle.
The fraction (e.g., 0.6) to be used for each
element was selected according to the proportion of each basic tissue in the composed
tissue. The values of the proportion used of each of the components were taken by
analyzing Fig. 5.2.
68
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Similarly^ the same approach was used to generate the bones, namely the ribs
and backbone. The available values of the dielectric properties of the cortical bone,
cancellous bone and bone marrow were used to generate the ribs and backbone. The
proportion of the components to simulate the bone was taken from [46].
Figure 5.3 shows the torso slice, with the organs within it. Chapter 6 reports the
results of the simulation of microwaves interaction with the torso, when a tumor is
inserted in the lungs. The tumor will be of different dimensions and located at different
positions.
There is not literature about the specific value of properties of a lung tumor. In
fact, most of the available literature and research focus in the breast cancer detection. In
this research area, researchers have demonstrated that the malignant breast tissues have
dielectric properties nearly identical to muscle, while normal breast tissues have
dielectric properties similar to those of fat.
fri this work the properties were defined using the physical constitution of the
tumor. Hence, it was also used a linear combination of the dielectric properties of the
following tissues: blood, blood vessel, and muscle. To adjust and verify the value of the
proportion of each of the corresponding components, the resulting dielectric properties
were compared to some o f the results reported at different frequencies [7] [47].
Figure 5.7 shows the dielectric properties of the tumor using the approach
mentioned above.
69
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
I - e - Relating pBrmittiwty
- 0 - Condyctivity IS/m)
Frequency (GHz)
Figure 5.7 Tumor dielectric properties for the frequency range 1-30 GHz.
As a next step in the process of computing the dimensions of the torso using
cells as distance units, the maximum required cell size at each tissue was found by
(5.1)
where
Sr i f ) = Relative permittivity in the tissue at frequency / .
Aqi f ) = Wavelength in free space at frequency / .
R = Number of cells per wavelength (10 were used for this thesis).
70
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
From the set of values of maximum cell size in each tissue, the smallest one was
selected, since is the one that meets the requirement of maximum cell size for every
tissue.
Once the value of Ax is obtained, and knowing the width and height of the
torso, the number of cells required to simulate the torso are computed.
In order to generate the corresponding computational array, an initial image of
the torso was created, whose dimensions are 794 by 543 pixels. According to the
required number of cells to build the torso, this picture has to be rescaled, it means that
either more cells (high-frequency analysis) or less cells (low-frequency analysis) may
be needed to represent the torso. For example, for a 1-GHz analysis the torso was
represented by 105 by 72 cells, while to represent the torso for a 10-GHz analysis, a 964
by 659 cells array was needed.
The percentage of either reduction or amplification of the picture is calculated.
Once this value is calculated, it is applied to the picture, using MS Paint application.
For instance, for a 10-GHz analysis, it is obtained from the software application that is
needed to increase the original picture by 21%.
This software application also generates a file containing the required
dimensions (in cells) of the complete FDTD space, which must be large enough to
contain the torso, as well as the measuring screen. The position of the total/scattered
field boundary is also contained. These values will be needed for the simulation later.
In addition to the dimensions file, this program also generates an input file,
containing some parameters for the simulation, such as the coordinates of both the
71
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
initial and final points for the screen, or measuring line. From this line the values of
are going to be recorded. Another parameters also contained in this file are: cell size,
frequency of analysis, number of iterations, etc. Due to the large amount of parameters
needed for each of the tests to be performed, it is very useful to generate automatically
this file, so it is ensured that the correct parameters are used in the simulation.
5.4 Converting the image to identifiers arrays.
As a final stage of the generation of the torso computational model, another
software application was developed [Appendix D]. This application converts the image
of the biological phantom into the tissues identifiers arrays.
Initially, the numerical values of the main colors available in MS Paint
application were obtained (RGB format).
Next, these colors are associated to each of the distinct tissues. For example,
free
space
is
associated
to
white
(R=255,0=255,B=255),
blood
to
red
(R=255 ,G=0,B=0), and so on.
Thereafter, the information of permittivity and conductivity of the distinct
tissues for a specific frequency were loaded.
A corresponding BMP file for each of the cases of study was generated. Hence,
according to the frequency of analysis, as well as the position and diameter of the
tumor, the corresponding file was analyzed.
72
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
m
103
230
an
4®
SOJ
603
Figure 5.8 Values of relative permittivity
) for each of the tissues for
frequency = 5 GHz.
An identification number, namely the identifier, was set to each of the tissues.
These identifiers will be used later in the simulation. Since integer numbers were used
to represent the tissues, the size of the file is not as big as it would be if instead of
saving identifiers numbers, the values of permittivity and conductivity for each of the
points of the grid were saved. For low-frequency analysis the file size may not be a
problem, but as the frequency of analysis is increased, the file size could represent an
inconvenience in terms of memory management. Besides, since a remote computer was
used for the simulation, the transfer time of the files was also an issue.
The main section of the program is the scanning of the image. The image is
read pixel by pixel. According to the color, the identification number of the tissue
associated to that color was saved in a file.
73
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
For illustrative purpose, an array with the permittivity values, which are
frequency dependent, was generated. Figure 5.8 shows the values of permittivity for
each of the tissues for a 5-GHz analysis.
WPI
t
I
w
u
Figure 5.9 Array containing the identifiers for a 5-GHz analysis.
Figure 5.9 shows the plot of the file containing the identification number for an
image 652 by 652, which is the image generated for a 5-GHz analysis.
After generating the identifiers file, another file containing the dielectric
properties for each o f the tissues for the frequency of interest was also generated. Table
5.1 shows the values of relative permittivity and conductivity for each of the tissues for
a 1-GHz analysis.
74
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table 5.1 Values of relative permittivity and conductivilcy for the tissues at IGHz.
Tissue
Relative permittivity
Conductivity (S/m)
Free space
1.000000
0.00000
Lungs
21.82500
0.47406
Blood
61.06500
1.58290
Bones
16.60784
0.26935
Heart wall
59.29000
1.28360
Interstitial fluid
68.87500
1.66730
Esophagus
64.79700
1.23160
Intercostal muscles
54.81100
0.97819
Thoracic wall
28.74150
0.50800
Tumor
56.51340
1.29111
A similar approach to generate the computational model is also used in [48] [49].
They carry out also an FDTD simulation, but to compute SAR and thermal elevation of
the human eye and head. Unlike our work, they just simulate at two frequencies.
75
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 6
SIMULATIONS, RESULTS, AND ANALYSES
6.1 In tro d u ctio n
The purpose of this research is to analyze the use of microwave for the detection
of a cancerous tissue in the lungs, by performing a numerical study. This study is
basically carried out by analyzing the effect of a malignant tissue in the scattered field,
as a response to an impinging electromagnetic plane wave.
The physical basis for tumor detection using microwave imaging is the
significant contrast in dielectric properties of background and malignant tissues. Figure
6.1 shows the values of these properties; namely the relative permittivity and
conductivity, o f lung [42] and a cancerous tissue [see section 5.3]. Although there is a
significant contrast in the entire range 1 - 30 GHz, this contrast is higher at low
frequencies.
In order to analyze a variety of cases, several configurations were set. This
way, it was possible to analyze the effect of a malignant tissue in the transmitted signal
through the torso, for different sizes and locations of the tumor. The simulations were
then performed using these different configurations.
The frequency of the incident
wave was also varied, in order to analyze a range of frequencies. This chapter presents
76
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the results obtained from the siimilations using the mentioned configurations, as well as
the analyses applied to them.
O
0
-f
Lung rslativs parmittlvity
Lung conductivity (S /m )
T um or raiatit® psrm ittivtty
T um or conducti's^y (S /m )
-
"*n
♦ ♦
i
U . . . ........... L......................
,
....
, + ■
i
i
25
30
F req u e n cy (GHz)
Figure 6.1 Dielectric properties of lung and tumor.
6.2 System configiiration
Before performing the simulations, a general configuration was set. Figure 6.2
shows the general configuration of the system used to analyze the effects of a tumor.
An EM plane wave strikes the torso model and the transmitted signal through the torso
is read at a screen located behind the torso.
It should be noted that the screen is located at the scattered field zone, in order
to be able to read the propagated signal through the torso model, without being affected
by the incident wave propagating around the torso. In Fig. 6.2 the model of a torso
without any malignant tissue is shown.
77
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Scattered/Total field boundary
Measuring line (Screen)
Absorbing
boundary
conditions
Torso
model
1
Incident plane wave
Figure 6.2 Diagram of the simulation o f an EM wave impinging
on the torso model.
For further simulations this initial model was substituted with the corresponding
model according to the size and position of the tumor. The dimensions of the simulated
torso are 38 cm. by 26 cm. The length of the screen is 46 cm.
In order to have several configurations, two positions and three sizes for the
tumor were established. The first position is located at the edge of the right lung, 14 cm
from the front surface of the model and 3.5 cm from the right side of the model; while
the second position is located at the inner edge of the left lung, close to the heart, also at
14 cm. from the front of the surface of the model, and 13 cm from the left side of the
model. The tumor was represented by a circle. Three sizes of the tumor were also
defined to use in the simulations. The selected diameters for the circle representing the
78
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
tumor were: 1, 3, and 5 cm. Figure 6.3 shows a 1-cm diameter tumor located at the
positions mentioned above.
m
I
1
ww^
Ia
I
r
Wr
b)
a)
Figure 6.3 Positions of the tumor for the study a) Position 1, at the outer edge of right
lung, and b) Position 2, at the inner edge of left lung.
6.3 Simulation approach
The method applied was the following: the signal on the screen was obtained
initially without any tumor present in the lung area. This signal is called reference or
background signal. Then, once a tumor is inserted, the simulation is performed. The
size and position and size were changed for each case of analysis. The reference signal
is subtracted from the resulting signal in each of the cases. The resulting signal from
this operation is called difference signal.
Hence, the main approach consisted of performing FDTD simulations for a) No
tumor present, to establish the background signal; b) a variable diameter of a tumor
having a fixed location; and c) change the location of the tumor and vary again the size.
79
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
To start the analyses, a suitable set of frequencies from the 1-30 GHz range was
set. The established set of frequencies to analyze the response was: 1, 5, 10, 14, 18 and
28 GHz.
Figure 6.4 shows the amplitude of the electrical field Ez at the screen for the
mentioned set of frequencies.
A plane wave impinging on the torso, at those
frequencies, was simulated. Analyzing the response, it can be seen that for the lowest
frequency of study, the response is very smooth, while for the highest frequency of
study; the response shows more the effect of the organs inside the torso due to the
decreasing of the wavelength.
^ 0.5
5 :g h z
D
0
-50
50
1
0.5
D
-
0
-500
0 .5
-S B
500
500
1
0.5
0
^
-5 0 0
0
500
0.5
-1000 -5 0 0
Cells along screen
500
1000
Cells along screen
Figure 6.4 Amplitude of the electrical field Ez read at the screen for the torso model
without tumor. The simulated frequencies are 1, 5, 10, 14, 18, and 28 GHz.
80
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Results were obtained for 1 and 5 GHz tests (shown in section 6.4), but for
higher frequencies the difference signal was so low, that it was not possible to find a
relationship between the difference signal and the size and location of the tumor.
Therefore, the research was concentrated for the range 1 - 6 GHz.
As initial step for the 1 - 6 GHz analysis, a simulation was performed to obtain
the response without a tumor present, that means, to obtain the background signal, also
referred to as reference signal. Figure 6.5 shows the transmitted signal through the
torso for each of the frequencies.
ty-
X
X
0.5
0^5
1 GHz
-50
___
yi 0.5 /
/
-50
1
\
i
0
-200 -100
1
i
IM. 0.5 /
/
\
0
200
Cells along screen
100
200
W '\
/
\
5 GHz
0
100
\ >,
L,
-200
i
V
4 GHz i
j/
I
50
0.5
100
-100
\
2 GHz
3 GHz
1
0
/
/
0
0.5
y
-100
50
1
\
N
m
6 GHz
-200
0
\
\
\
200
Cells along screen
Figure 6.5 Amplitude of the electrical field Ez read at the screen for the torso
model without tumor. The simulated frequencies are 1, 2, 3, 4, 5,
and 6 GHz.
81
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
6.4 Steady State Analysis
The study was divided in two approaches: a) Steady state response, and b) Time
series analysis. For the steady state response it was found the difference between the
reference signal and the resulting signal after inserting a tumor, once the system had
reached the steady state. With respect to the time series analysis, the response of the
system was analyzed paying attention to the transient behavior of the difference signal.
The first analysis performed was the steady state response of the system.
According to the defined positions and dimensions for the tumor, several cases were
analyzed. Table 6.1 lists these cases. The signal at the screen is saved, and thereafter,
the background signal is subtracted from them.
Case
1
2
3
4
5
6
Table 6.1 Description of cases to analyze
Diameter of
Position of tumor
tumor (cm)
1
Outer edge of right lung
Outer edge of right lung
3
Outer edge of right lung
5
1
Inner edge of left lung
3
Inner edge of left lung
Inner edge of left lung
5
C A SE l
Diameter of tumor: 1 cm
Position: Outer edge of right lung
82
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 6.6 shows the configuration of the first case of study. An incident plane
wave striking on the torso slice, with a 1-cm tumor located at the outer edge of the right
lung.
Measuring line (Screen)
/
Incident
plane wave
Figure 6.6 Diagram of case 1. EM plane wave impinging on the torso slice
with a 1-cm diameter tumor located at outer edge of right lung.
The simulation was performed for all the frequencies of interest, namely, 1 - 6
GHz. After reaching the steady state, the amplitude of the signal in the screen was
computed. Figure 6.7 shows the corresponding values of the difference signal (AEz).
As a reminder, the shown signal is the difference between the read signal with a tumor
present, and the signal without a tumor on the model.
In Fig. 6.7 it can be seen that for low frequencies the signal difference is
relatively large, but the position of the signal peak with respect to the tumor axis is not
as accurate when compared to high frequencies. However, for high frequencies the
level of the difference signal is low. Table 6.2 shows the values of the signal peaks for
each of the frequencies.
83
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
r ' \
........
A ........
A
. \ ---
'
ol...
"V ------
......
. . .
-----
;■ ” r
xIO"*
r
~JT TT
T
T
Figure 6.7 Difference signal for case 1. EM plane wave impinging on the torso
slice with a 1-cm diameter tumor located at outer edge of right lung.
Table 6.2 Values of maximum peak of the difference signal for case
Frequency
Peak AEz value
(GHz)
(V/m)
(dB)
1
5.5 X 10'^
-35.29
2
2.6 X 10^
-35.85
3
4.25 X 10’^
-43.72
4
1.1 X 10'^
-49.58
5
2.0 X 10"®
-57.00
6
0.4x10-®
-64.00
84
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CASE 2
Diameter o f tumor: 3 cm.
Position: Outer edge of right lung.
Measuring line (Screen)
/
Incident
plane wave
1
1
1
1
1
1
1
Figure 6.8 Diagram of case 2. EM plane wave impinging on the torso slice
with a 3-cm diameter tumor located at outer edge of right lung.
The configuration of the second case of study is shown in figure 6 .8 . An
incident plane wave striking on the torso slice, with a 3-cm diameter tumor located at
the outer edge of the right lung.
Similarly to the first case, the simulation was also performed for each of the
frequencies of interest.
Figure 6.9 shows the results of the simulation. As in case 1, the value of the
peak of the signal is greater for low frequencies, but the position of this peak is more
related to the location o f the tumor for high frequencies. The values of the difference
signal peaks are listed in Table 6.3.
85
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
j
;
;
A
/
\
/
^ .....
; / 1
\
: /
!
1
; _/__:
-S3
i
■
\
tm
-40
Tum<Jax«.
-20
—
0
20
Cetis atong ^cresr.
43
-100 t
-50
Tumor asis
60
1
\
0
CcHs along scfaen
1
1
SO
100
1
\
\
\
/
/
\
/
/
-150 t -100
Tumor sxSs
5o
0
cstis along s<
4c a r
i
•
SO
Tm
•253
IKJ
-200 -150 -100
TumiJastB
-300
T-20O
i
-50
0
SO
CsKs slang screan
-1O0
0
100
100
IK)
200
200
230
300
Celts along
Figure 6.9 Difference signal for case 2. EM plane wave impinging on the torso
slice with a 3-cm diameter tumor located at outer edge of right lung.
Table 6.3 Values of maximum peak of the difference signal for case 2.
Frequency
Peak AEz value
(GHz)
(V/m)
(dB)
1
5.5 X 10'^
-32.60
2
2.0 X 10"^
-37.00
3
7.5 X 10'^
-41.25
4
0.9 X 10'’
-50.45
5
2.5 X 10-®
-56.00
6
0.5 X 10'®
-63.00
86
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CASES
Diameter of tumor: 5 cm.
Position: Outer edge of right lung.
Measuring line (Screen)
/
I
/
t Tt t t TT
Incident
plane wave
Figure 6.10 Diagram of case 3. EM plane wave impinging on the torso slice
with a 5-cm diameter tumor located at outer edge of right lung.
Figure 6.10 shows the diagram of the third case of study. A tumor with a 5-cm
diameter is located at the outer edge of the right lung. The incident plane wave goes
through the torso model and the scattered signal is read at the screen.
The difference signals for each of the frequencies are shown in Fig. 6.11. The
value of the signal peak is greater, as expected due to the increment in the tumor size.
Table 6.4 shows the peak signal values.
87
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
/
/
!
/
\
\
!
;
1
\
\
-so
-20
TurrsJ a>03
1
'
!
:
2GH j
-100
0
20
“‘“ ’3 screfiti
1
;
X
1GH*
---^
x "' i
—
;
\
/
i /
i/
:
:
T -50
Tumor axis
-
0
Calls atongsoraon
50
1
1
\
...1
1
/
i
.
i
J . ...
XT--.
ui“ °
\
i
/
—
•1S0
-130
Turno'r ante
-SO
....
. . . .
3GH*
\
- / -
0
atwig si
0
-ISO -100
t
Turner axis
!
-SO
0
50
Cells along screen
I
!
100
!
ISO
--
200
1
...
...
I
3
-2t»
sa fe
-1K)
0
1CO
Celte along soraon
200
300
Figure 6.11 Difference signal for case 3. EM plane wave impinging on the torso
slice with a 5-cm diameter tumor located at an outer edge of right lung.
Table 6.4 Values of maximum peak of the difference signal for case 3
Frequency
Peak AEg value
(GHz)
(V/m)
(dB)
1
8.0 X 10-^
-31.00
2
3.0x10-^
-35.23
3
8.0 X lO"’
-41.00
4
1.4 X 10'^
-48.54
5
3.5 X 10-®
-54.56
6
0.3 X 10'®
-65.23
88
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CASE 4
Diameter of tumor: 1 cm.
Position: Inner edge of left lung.
Measuring line (Screen)
/
Incident
plane wave
t t t t t Tt
Figure 6.12 Diagram o f case 4. EM plane wave impinging on the torso slice
with a 1-cm diameter tumor located at inner edge of left lung.
The configuration for the fourth case is shown in Fig. 6.12. Contrary to the
three previous cases, the tumor is now located at the inner edge of the left lung, close to
the heart. The size of the malignant tissue is 1 cm.
Similarly to the previous cases, the simulation was performed for each
frequency of interest.
Figure 6.13 shows the results of the simulations.
Some special results are
obtained for this case. Contrary to the first three cases, in which for each frequency
there was a predominant peak, in this case there are two frequencies, in which the
difference signal does not have a predominant peak (4 and 6 GHz), while in another (5
GHz) the signal has an irregular pattern.
89
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
;
-
;
/
r-'-
%um
-so
-2Q
0
Cdle
20
40
sciesn
__rrrrr,-
j ! ; 1
; ! ; 1 1
/ i \
\ ;
i
\-V- // ' - \ - ./---A__
/ i
......
/
-1Sfl
-100
-50
0
SO
Calls siona sMSsn
•
ICO
•2S0
1S3
-200
-ISO
-ICa
-50
0
50
Cells al«ig screei^
100
ISO
t
^
200
------i.......
r>Hr-3 ?
----- y-f-- ] JJV....
-303
-200
-100
-300
-200
-10Q
0
100
Cells along e c r e ^
^
200
SCO
Figure 6.13 Difference signal for case 4. EM plane wave impinging on the torso
slice with a 1-cm diameter tumor located at inner edge of left lung.
A special result was expected due to the proximity of the tumor to hearth.
These two tissues have similar dielectric properties. Besides, the penetration depth
decreases as the frequency increases, and the random scattering within the torso
increases as the wavelength is decreased. Hence at high frequencies it is not possible to
detect a predominant peak.
90
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Table 6.5 shows the values of the signal peaks. When compared with case 1, in
which the tumor has also a diameter of 1 centimeter, the values are lower. Thus, the
effect of the heart is evident, due to the similarities of its electrical properties with those
of the cancerous tissue.
Table 6.5 Values of maximum peak of the difference signal for case 4
Frequency
Pealc AEz value
(GHz)
(V/m)
(dB)
1
5.5 X 10"*
-32.60
2
0.2 X 10'^
-47.00
3
2.0 X 10’®
-57.00
4
No predominant peak
5
1.3 X 10'*
-58.86
6
No predominant peak
—
CASES
Diameter of tumor: 3 cm.
Position: Inner edge of left lung.
Figure 6.14 shows the configuration of the fifth case of study. The incident
plane wave strikes on the torso slice, with a 3-cm diameter tumor located at the inner
edge of the left lung.
91
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Measuring line (Screen)
/
Incident
plane wave
t Tt t t Tt
Figure 6.14 Diagram of case 5. EM plane wave impinging on the torso slice
with a 3-cm diameter tumor located at inner edge of left lung.
Similarly to case 4, only for the low frequencies (1, 2, and 3 GHz) there is
present a predominant peak, while for the rest of frequencies there is not such a peak.
Figure 6.15 illustrates that when the incident wave has a frequency of 4 and 6
GHz, there is not a predominant signal peak, while for 5 GHz the pattern of the
difference signal is irregular. The values of the difference signal peak are listed in
Table 6.6 .
92
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
: /'
; ;
/<';......f!
i «
1
;
!
:
;
;
1
!
i
i
!
I
1
....I..;.,....i .. J j.i.,..;,
.....
.....
...
1 :
:
\
.....
.....
3. . .
4 G Ife
....; .
.....
Calls along S!
J 0
-300
CellB along:
-2C0
-100
0
' 100 T 200
300
Celfe Btong sereanT*"™^ ®*®
Figure 6.15 Difference signal for case 5. EM plane wave impinging on the torso
slice with a 3-cm diameter tumor located at inner edge of left lung.
Table 6.6 Values of maximum peak of the difference signal for case 5
Frequency
Peak AE^ value
(GHz)
(V/m)
(dB)
1
1.55 X 10'^
-28.23
2
0.8 X 10'^
-41.00
3
4.5
-53.47
X IQ -*
4
No predominant peak
5
0.6 X 10'*
-62.22
6
No predominant peak
—
93
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CASE 6
Diameter of tumor: 5 cm.
Position: Inner edge of left lung.
Measuring line (Screen)
/
Incident
t t TTt t t
Figure 6.16 Diagram of case 6. EM plane wave impinging on the torso slice
with a 5-cm diameter tumor located at inner edge of left lung.
The configuration of case 6 is shown in figure 6.16. An incident plane wave
striking on the torso slice, with a 5-cm diameter tumor located at the inner edge of the
left lung.
To find the response in the screen, the simulations were performed for each of
the frequencies in the range 1- 6 GHz. Results (Fig. 6.17) show that the steady state
signal has a similar behavior as in cases 4 and 5. That is, the signal either does not
show a predominant peak or does not have a regular pattern for certain frequencies.
Table 6.7 shows the values of the signal peaks, when present.
94
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
7
:
/
/
/
/ \
I/—-:
=■
/
/
.... \
fii -
•00
-40
-20
-
/
2GHz
,
0
20
40
C9!!s aiana scraan Tunwsxta
Is alor^ screen
;
;
;
;
:
;
IS8W
-ISO
-lO )
-«0
0
so
105
CeUs along scfeen ju n v J acts
150
■2S0
y
i
i
i
<Gfb
;
;
i
-200
-130
-100
-50
0
50
Celte along screan
ICO
t
150
200
250
(A,
• ¥ \i
^
•3C0
-200
-100
^ 0
~
1OT t
2£0
C a itsaiw s screen Tumor asb
200
^
------- ------------
-sra
__
0
io o
1 2M
is along screan Tumor a
---------
am
Figure 6.17 Difference signal for case 6 . EM plane wave impinging on the torso
slice with a 5-cm diameter tumor located at inner edge of left lung.
Table 6.7 Values of maximum peak of the difference signal for case 6 .
Frequency
Peak AEj value
(GHz)
(V/m)
CdB)
1
2.00 X 10'^
-27.00
2
1.0 X 10"*
-40.00
3
1.8 X 10'^
-47.45
4
1.0 X 10'®
-60.00
5
1.5 X 10'*
-58.24
6
No predominant peak
—
95
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
«■,
/
\
W S
/ \ i
i
-.1
\
\
i\ ;
; \ ;/
i '-T
af »|J
'v / i
/ : \ \ ^;
i
i
;.....
\\
\
;
;
si
i"''''
:
i
:i 1
^
'
' \
:
;
/"K
/ ; h-T
/ ‘-G— ....:---
Figure 6.18 Difference signal (AEz) for 3-GHz steady state analysis. The
tumor is located on outer edge of right lung.
Several cases were analyzed. In each of them either the size or the position of
the tumor changed.
Each case was simulated using a different frequency for the
incident wave. Figure 6.18 shows the difference signal when an incident wave of 3
GHz strikes the torso. The tumor is located at the outer edge of the right lung.
I
^r.
1
\
/
/
/
/
..A ...:../
/N
‘vy . ;'L. (
kk
h
\ /
/ ^
;\
; ^
I
1
\
;
!
sc
l
V i/
; V
]
\
\
^ I''^-4 :
:
!
100 160
Figure 6.19 Difference signal (AEz) for 3-GHz steady state analysis. The
tumor is located on inner edge of left lung.
96
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 6.19 shows the results wheu the position of the tumor is changed to the
inner edge of the left lung. In every case there is a predominant peak, and its position is
closely related to the location of the malignant tissue. Hence, one of the conclusions of
the study is having found, among the analyzed frequencies, what seems to be the
appropriate frequency for lung tumor detection. The criterion to select it was based on
two parameters.
The first parameter is the value of the predominant peak of the
response, if present. The second parameter is the relative location of this peak with
respect to the location of the tumor. Using these criteria, and based on the results, 3
GHz appears to be the appropriate frequency for lung cancerous tissue detection, since
it is the one, which presents a predominant peak at every case. The position of the
signal peak can also be related to the position of the malignant tissue.
6.5 Time series analysis
One of the advantages of using a time domain method to simulate the
interaction of electromagnetic waves with a target is the fact that the wave propagating
can be visualized at every time step. This way, it is possible to find and analyze any
transient behavior of the wave; that means, any phenomenon before reaching the steady
state. In this work an approach called “time series analysis” is presented.
This analysis is performed on a plot of the transmitted field through the torso,
arriving to the screen at each time step.
The shown values are the result of the
subtraction of the background signal from the signal obtained when a tumor is present.
Figure 6.20 illustrates how this concept is accomplished.
97
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Screen
m
Plane
Wave
iw
Si
■
ti
ty_
Figure 6.20 Diagram of the time series analysis. The values at the screen
are shown at every time step.
This analysis was performed for each of the cases, at each of the frequencies of
interest; however, only a few representative cases are illustrated here with the purpose
of showing the information that can be obtained by analyzing the transient behavior of
the scattered signal.
Figure 6.21 Time series of scattered AEz at the screen for a 1-GHz incident
wave. A 1-cm tumor is located at outer edge of right lung.
98
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 6.21 shows the time series progression of the difference signal at the
screen, when a 1-GHz plane w'ave strikes the torso model. A 1-cm tumor is located on
the outer edge of the right lung. The cells representing the screen are located on the
vertical axis.
The larger difference occurs between 3 and 4.5 ns, after starting the
simulation. This suggests that analyzing the signal, at this time, could offer the best
time to detect the cancerous tissue. This can be experimentally implemented by using
Time Domain Reflectometry (TOR).
Similarly, Fig. 6.22 shows the time series of the difference signal at the screen
when a 1-GHz plane wave strikes the torso model. In this case a larger tumor, 5-cm
diameter, is located on the outer edge of the right lung. The larger difference takes
place between 3 and 5 ns, after starting the simulation.
ff
S
10
12
14
7ifrw (ps)
Figure 6.22 Time series of scattered AEz at the screen for a 1-GHz incident
wave. A 5-cm tumor is located at outer edge of right lung.
99
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 6.23 shows the time series of the difference signal at the screen when a
l"GHz plane wave strikes the torso m odel In this case, a 1-cm tumor is located at the
inner edge o f the left lung. The difference in the early time is not as evident as in the
previous cases, when the tumor was located on the outer edge of the left lung. This is
due to the proximity of the tumor to the heart. This behavior was expected, since the
dielectric properties of the heart and the cancerous tissue are similar.
Besides, the
relative big size of the heart, with respect to the tumor, and its high conductivity makes
the heart act as a relatively good scatterer, which fades the effect due to the tumor.
o
80
a
M
12
14
Time ins)
Figure 6.23 Time series of scattered AEz at the screen for a 1-GHz incident
wave. A 1-cm tumor is located at inner edge of left lung.
In the case presented in Fig. 6.24, the incident field has a frequency of 1 GHz.
The position of the malignant tissue is also on the left lung, however its size is 5-cm.
Due to this increment in size, it is possible to detect a difference in the signal, which is
evident at 4 ns.
100
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
.•-J'
2
4
Q
S
10
Time(rw^
12
14
16
ia
20
Figure 6.24 Time series of scattered AEz at the screen for a 1-GHz incident
wave. A 5-cm tumor is located at inner edge of left lung.
Figure 6.24 also shows that after approximately 8 ns the system reaches its
steady state.
Figure 6.25 Time series o f scattered AEz at the screen for a 3-GHz incident
wave. A 1-cm tumor is located at outer edge of right lung.
In the case presented in Fig. 6.25, the tumor is located at the outer edge of the
right lung. The size of the tumor is 1 cm. The incident field has a frequency of 3 GHz.
101
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The significant difference is present between 2 and 3 ns. After approximately 4 ns the
system reaches its steady state.
Tt^r»r ->
Figure 6.26 Time series of scattered AEz at the screen for a 3-GHz incident
wave. A 1-cm tumor is located at outer edge of right lung.
Figure 6.26 shows the difference signal, for the same case, between 18 and 24
ns. It is evident that the system has already reached the steady state.
MS—
I
a»s
Ttmefra)
Figure 6.27 Time series of scattered AEz at the screen for a 3-GHz incident
wave. A 1-cm tumor is located at inner edge of left lung.
102
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The period of analysis for the case shown in figure 6.27 would be between 3 and
4 ns. The size of the simulated tumor is 1 cm. A 3-GHz plane wave impinges on the
torso. The tumor is also located at the inner edge of the left lung, meaning that it is
close to the heart.
50
ISTo?
150
r
je
§ 330
SSO
iB i
s®Mi
Figure 6.28 Time series of scattered AEz at the screen for a 4-GHz incident
wave. A 1-cm tumor is located at outer edge of right lung.
For frequency 4 GHz, even in the early time, the tumor is only detectable when
it is located away from the heart (Fig. 6.28), whereas when it is positioned close to the
heart (Fig. 6.29), the tumor is not detectable.
103
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
0^
1
1.5
2
2 .5
3
3.5
Tima (n8)
Figure 6.29 Time series of scattered AEz at the screen for a 4-GHz incident
wave. A 1-cm tumor is located at inner edge of left lung.
Hence, when the tumor is located at the inner edge of the left lung, the
simulation shows that it is not detectable, that means, there is not an evident difference
even at early time. This is valid for 4-, 5-, and 6-GHz analysis.
X10
BHHHI
•2
1 00
i1.5
200
P
l«o
i
500
■
9 H H II
700
"is"
S
35
Time (nsi
Figure 6.30 Time series of scattered AEz at the screen for a 6-GHz incident
wave. A 1-cm tumor is located at outer edge of right lung.
104
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 6.30 shows that for a 6-GHz wave, a 1-cm tumor is visible, although
slightly, only when the tumor was located away from the heart.
The cases illustrated in this section show that the analysis of the time series
could be useful to find the optimum time, during the incidence of a EM plane wave, to
detect a malignant tissue.
As in the steady-state analysis, it is difficult to detect a
difference when the tumor is located in the proximity of the heart.
6.6 Conclusions
One of the conclusions after analyzing the results from the simulations is having
found the appropriate frequency for detection of a malignant tissue in the lungs. The
criterion to select the frequency was based on two parameters. The first parameter was
the value o f the predominant peak of the response, if present. The relative location of
this peak, with respect to the position of the malignant tissue, was the second parameter.
Using these criteria, and based on the results shown in this chapter, 3 GHz appears to be
the optimum frequency for lung cancer detection, since it is the one, which presents a
predominant peak at every case. Moreover, the position of the signal peak is closely
related to the location of the malignant tissue.
In addition, the time series analysis can be useful to find the optimum time to
detect a malignant tissue. Results also shown that either for steady-state analysis or
time series analysis, it is difficult to locate a difference when the tumor was located in
the nearness of the heart. This issue is due to the similarities of the electrical properties
of the heart and those of the cancerous tissue, as well as the difference in size.
105
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
CHAPTER 7
CONCLUSIONS AND FUTURE W ORK
7.1 Conclusions.
Microwave systems for cancer detection are being researched. Enough promise
is there to be optimistic that a viable system will become available. It is not clear if
such a system or systems will replace x-ray chest detection as a screening tool for lung
cancer. But a combined use of these systems can improve detection and limit false
positive findings. In addition, microwave technology is nonionizing and noninvasive.
In the present work, the results of a simulation of a plane wave impinging on a
torso slice are presented.
Simulations were conducted using several frequencies for
several sizes and locations of a cancerous tissue in the lung area. The dimension and
the position of tumor were varied in order to find the most suitable frequency in the
range 1 - 3 0 GHz for the lung cancer detection. It was found that there were significant
background scattering from lung and other tissues above 6 GHz.
focused on the range of 1 - 6 GHz.
Hence, the study
For low frequencies in this range, the signal
difference between the background signal and the signal obtained by inserting the tumor
was relatively large, but the information of the position of the tumor is not as accurate
when compared to high frequencies. The drawback of the high frequencies analysis
was the low level of the difference signal due to an increase in background scattering.
106
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The reason of these results is that the spatial resolution of a small tumor
improves at shorter wavelengths. However the depth of wave penetration, and therefore
the amplitude of the scattered response of the tumor, drops at shorter wavelengths.
Considering this trade off between level of the signal and location of the peak
signal with respect to the position of the tumor, the frequency 3 GHz appears to be the
optimum frequency for detection of lung cancer.
Difficulties were found to detect a tumor near to heart due to similarities in the
electrical properties, and the relative big size of heart with respect to small tumors. The
random scattering of the torso tissues and the skin-interface reflection frirther obscure
the scattered response of a small tumor.
A time series analysis is also presented as a method to find the best time to
measure the highest difference in the signal. The difference signal in several time steps
were presented in this report and the optimum time to detect the presence of a tumor
was discernible.
7 .2
Future w ork.
The work presented here is an initial study of the use of microwave frequencies
for detection of cancerous tissue in the lung. This study showed a potential for lung
cancer detection using a microwave system.
Recommendation for future studies include (1) applying a more advanced data
processing to the measured signal; (2) performing the simulation using different system
configurations, for example, rotating either the torso model or the angle of incidence of
107
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
the impinging wave. By doing this rotation, there should be enough data to generate the
image; (3) investigating the use of not only a plane wave as incident field, but also other
type of sources in order to be able to focus the incident field in a specific area; (4)
developing and utilizing 3-dimensioal and more realistic organ/tissue computational
models; and (5) developing and conducting experimental studies using tissue phantoms.
108
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX A
APPROXIMATIONS FOR DERIVATIVES
109
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Approximation for derivatives.
The following description of approximation for derivatives is mainly based on:
R.C. Booton Jr., “Computational Methods for Electromagnetics and Microwaves”,
Wiley Series in Microwave and Optical Engineering, 1992.
To start the description, consider the Taylor expansion
(A.1)
f ( x + h) = f ( x ) + hf'(x) + y / " W + ■•
from which
(A.2)
2
h
/ 'W '
f{x +h)-f{x)
(A.3)
Equation (A.3) is referred to as Forward difference formula. Using the same
procedure it can be obtained
jy y /W
z /k l^
(A.4)
which is referred to as Backward difference formula. These are an approximation to the
derivative /'( x ) , but as h approaches 0, the error approaches 0 only as the first power of
h and thus is not a desirable approximation.
The left side of (A.2) is a better
approximation to the derivative at the midpoint. To see this, consider the expansions
f{x +h ) = f
h
X +
V
/
/2^
+ - / ' X 4" +
2j 2 V
2,
8
r
rx + —
h) + —
h
K
2)
48
/
r
X +
\
2)
(A.5)
and
/W =/
8
48'
110
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(A.6)
Subtraction of these expansions yields
f { x + h) ~ f { x ) _
^
h'^
x + — + - !_ /* x + ~- “h ‘
2y
V
24
V
(A.7)
2y
Note that the error decreases as the second power of h and hence this is a more
desirable approximation.
f{x + h)-f{x)
f\x +2)
(A.8)
h
This is a good approximation for the first derivative between nodal points, but
usually approximations at the nodal points are needed. To find such an approximation,
combine the expansion
f { x - h ) = f { x ) - hf'{x)^
2
O
+■
(A.9)
with (A .l) to get
(A.10)
Ar
\ f h = — , these expression can also be written as
/
Ax
x+- /
Ax
-
Ax
• ••
(A.11]
Hence,
Ax
/ ' W * / X+ '
- / I
X-
Ax
Ax
(A.12)
This approximation is referred to as central difference formula, which in general
it is preferred because o f the greater accuracy.
I ll
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX B
GRAPHS OF D IELECTRIC PROPERTIES OF BIOLOGICAL TISSUES FOR
FREQUENCY RANGE 1 - 3 0 GHZ
112
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
- S - RBlativa permittfnty
- t - Canductivnly (S/m)
FreqiiSficy (u H 2j
Figure B.l Lung dielectric properties for the frequency range 1-30 GHz.
Relative permittivity
-■€- CDnductiwty (S/m)
Figure B.2 Blood dielectric properties for the frequency range 1-30 GHz.
;
;
1
S
to
—^
1
15
20
Frequsricy (GHz)
Relative parmrltivHy
Condudwity (S/mj
i
25
Figure B.3 Bones dielectric properties for the frequency range 1-30 GHz.
113
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
F raqusncy (GHz)
Figure B.4 Heart wall dielectric properties for the frequency range 1-30 GHz.
Figure B.5 Interstitial fluid dielectric properties for the frequency range 1-30 GHz.
Figure B .6 Esophagus dielectric properties for the frequency range 1-30 GHz.
114
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure B.7 Intercostal muscles dielectric properties for the frequency range 1-30 GHz.
Figure B .8 Thoracic wall dielectric properties for the frequency range 1-30 GHz.
Figure B.9 Tumor dielectric properties for the frequency range 1-30 GHz.
115
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX C
MATLAB PROGRAM “TORSO DIMENSIONS IN CELLS”
116
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
%Torso dimensions in cells.
%This program find the minimum cell size to be used to build the
%striM:;ture according to the fiequency o f analysis.
%Since the higher the frequency, the lower the wavelength, there
%is going to be needed more cells to represent the same width, and
%height.
clear all;
^rint^l,'Start analysisVi);
fr'eq = input(Enter frequency to analyze(G H z):');
fprintf(l,'Frequency = %d GHzVi',fi:eq);
min_dx = mm(dx_tissues(freq,:));
fprintfr^ 1,T he minimum value at freq = %d GHz is% 1 0 .9 f
m'\n',freq,miB_dx);
cic;
%Dimensions o f torso in meters.
Q3rintf(l,'Dimensions o f torso in ceUs.W);
4)rintf{i,'Dimensions o f torso in meters.\n*);
^printf(1
width_cells = widtb/min_dx;
heightjcells = height/mm_dx;
rounded_width_cells = round(width_ce!ls);
rounded_height_celk = round(heig)it_ceiis);
,'Width: %d cells Height: %d
cells\n',rounded_width_cells,rounded_height_celis);
w id th - 0 .3 8 ;
height = 0.26;
^rintf^l,'Width = % 4.2f mts.
Height - % 4.2f
rats.\n',width,height);
GHz = le9; C = 3e8;
%The dimensions o f the original image w ere clianged
%to have the same ratio as 38/26.
%Previous width == 810 -> 794
%Height w as unchanged
%The dimensions o f the original image.
% (tQrso_flmHf_newwidth) are:
LUNGS = 1;
BLOOD = 2;
BONES = 3;
HEART__WALL = 4;
INTERSTTnAL^FLUID = 5;
ESOPHAGUS = 6;
INTERCOSTAL_MUSCLES = 7;
THORACIC_WALL = 8;
[lungs_freq_range,biood_freq__range,bones_freqj'ange,heart_waii_fieq_range,...
interstitial_fluid_freq_range,esophagus__freq_range,...
intercostal__muscles_freq_range,thoracic_wall_fieq_range] = ...
tissues_properties_freq_range;
^rintf(l
^rmtfr^l,Torso dimensions in the original BMP picture.W);
fprintf( |
width_orig_celis = 794; beight_orig_cells = 543;
^rinti(l,'Original Width: %d cells
Original Height: %d
cells\n',width_orig_cells,height_orig_celis);
fo ri= l:3 0
trequency_range(i) = i^GHz;
iarabdaO(i) = C/frequency_range(i);
end
EPSR =1;
%These dimensions correspond to the torso only, and not to the
%BMF file
SIGMA = 2;
%width_orig_cells * factor = rounded_width_cells;
%We need to find the "factor" value to use it
%in M S paint to rescale the original picture
%eps r
= Relative permttivity in the tissue at frequency freq
(GHz)
%larri>da = Wavelength in the tissue
%dx_tissues = Minimum cell size in the tissue vs. frequency
^rm tfl^l.R escalii^ factors.\n');
for freq = 1:30
eps_r = lungs_freq_range(freq,EPSR);
lambda(fieq) = lambdaO(freq)/sqrt(eps_r);
dx__tissues(freq,LUNGS) = lambda(freq)/IO;
eps_r = blood_freq_range(fi^,EPSR);
lambda{freq) = lambdaO(freq)/sqrt(eps_r);
dx_tissues(freq,BLOOD) = lambda(freq)/IO;
eps_r = bones_freq_mnge(freq,EPSR);
lambda(freq) = lambdaO(freq)/sqrt(eps_r);
dx_tissues(freq,BONES) = !ambda(freq)/10;
eps_r = heart_waU_freq_range(freq,EPSR);
lambda(freq) = lambdaO(freq)/sqrt(eps__r);
ds_tissues(freq,HEART_WALL) = to b d a (fr e q )/10;
eps_r = interstitial_fluid_fi-eq_range(freq,EPSR);
larabda(freq) = lambdaO(freq)/sqrt{eps_r);
dx_t!ssues(freq,lNTERSTJTIAL_FLUID) = lambda(freq)/10;
eps_r = esophagus_freq_range(freq,BPSR);
larabda(freq) = iarabdaO(freq)/sqrt(eps_r);
dx_tissues(freq,ESOPHAGUS) = lambda(freq)/10;
eps_r = intercostal_musc!es_freq_range(fr:eq,EPSR);
lambda(freq) = lambdaO(freq)/sqrt(eps_r);
dx__tissues(freq,INTERCOSTAL_MUSCLES) = lambda(freq)/10;
eps_r = thoracic__wali_freq_range(freq,EPSR);
lambda(freq) = lambdaO(freq)/sqrt(eps r);
dxJissues(freq,'m ORACIC_W ALL) = iambda(freq)/iO;
end
fectorl = rounded_width_cells / width_orig_cells;
factor2 = rounded_height_celis / height_orig_celjs;
frctori j j e r c = factorl 100;
factor2_perc = factor2 100;
roiinded_factorl_perc = round(factorl jserc);
rounded_factor2_perc = roimd(frctor2j>erc);
^rintfr^l,'Factor for Width and Height in percentage: % 2.2f
\n',rounded_factorl_perc);
%^rintfl[i,'Factor for Height in percentage: % 2.2f
\n',rounded_factor2_perc);
^rintfr[l
^rm t^ 1.Distances.\n');
%Distance in meters
%This value must not be less than 0.07 mts to ensure
%the torso is inside the total field space.
dist__torso__scat_field = 0.07;
dist_scat_fieM__ml = 0 .0 !;
d istjn!_p m l = 0.02;
pm l_cells = 1 5 ;
lprlntf(l,'Distance torso->scattered field = %f
mts.'‘>n’,dist_torso_scat_fieid);
^rintfr^ls'Distance scattered field->measuring line = %f
mts.\n',dist_scat_fleld_ml);
^rintfi;i,'Distance measuring line->PM L = % f mts.\n',dist_ml_pnil);
%Showing minimum cell size (dx) per wavelength.
dist torso scat field cells = dist torso scat field/min dx;
117
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
yp = rounded_xf;
xp = rounded_yf;
^rintfCi.'Vaiues in Paint coordioatesVi’);
fprint^ljxp = %d yp = ®/od\a'.xp,yp);
d!st__scat_field_jnl__cells = dist_scat__field_ml/ram_cix;
dist_mi_pml_cells == dist_ml_pirii/miii_dx;
^rintf(i,'Distance torso->scattered Held = %d
ceils.'«i’,rotind(dist_torso_scat_fieid_cells));
§)rinti( I .D istance scattered fieid->measuring line = %d
ceils. W ,round(dist_scat_fieldjnl_cells));
^rintS^ I,'Distance measuring line->PML = %d
ceils. W,round(dist_ml_pml_ce!ls));
^riiit^l^'FML = %d cells\rf,pml_cells);
%Ca!culate how many time steps for time
dt = min_dx / { 2 ^ C);
steps - input(Enter iterations:');
totai_time = steps * dt;
fprintf{l,'%d Iterations = % f nanosecs.\n',steps,tota]__time*le9);
time = inpu^'Enter time to simulate:');
steps = time / dt;
^rintf(l,'Required iterations = %d steps'®’,round(steps)};
^rintf(l,'FD TD dimeBsions\n');
JE = height_celis + (2*dist_torso_scat__field__ceils) 4-...
(2'*dist_scat_fieid_ml_cells) + (2®dist_ml_pml_cells) -f ...
{2*pml_cells);
IE = JE;
Generate_dimensioi^_h(rounded_IE,rounded_JE,IC,JC,roimded_IA,
rounded_JA);
rounded_IE = round(IE);
rounded^JE == round(JE);
fprintfCI.IE = %d
to = 25.0;
steps ~ inputCEnter iterations for input file :');
movie = input(Do you want to save E z for M ovie?(yes=l ,no=0): ’);
movie_period = 0; %De&ult value
if{ m o v ie ~ l)
raovie_period = inputCMovie p eriod:");
JE = %d\n’,roundedJE,rounded_JE);
%Checking whether IE and JE are even,
res = Toimded_IE/2;
ei^
if (res= fioor(res))
^rintf^l.'It is even!\n’,pml_celis);
else
^rintf(I,'It is oddi‘\n',pml_celis);
rounded_IE = rounded_IE + 1;
roimded_JE = rounded_JE + 1;
end
fprintt( I ,'IE = %d
spread — 8.0;
axis = 1; %Aiong X-axis
Generate_input_txt(min_dx,freq,tO,spread,...
pm l_cells,m ovie,mo vie_period,steps,...
axis,INITX,rounded_INITY,EKDX,rounded_ENDY);
JE = %d\n',rounded_IE,rounded_lE);
%Function to load the properties o f the tissues in the range 1-30
%GHz.
IC = roimded_IE/2;
JC = rounded JE/2;
ip rint^ l/IC = “%d JC =
mnctioa
[lun^_fr,blood_fi‘,bones_fr,heait_wall_fi-,mteretitial_flmd_fi',...
ssophagus_ff,intercDStal__miiscles_fr,thoracic_wall_fr,tumor__fr]
Total/Scattered field boimdary'o');
JA = pml_cells + diEt_ml_pmi__celis + dist_scat_field_ml__celis;
rounded_JA = roimd(JA);
roimded_IA = rounded_JA;
:^rintl{ 1,1A = %d JA = %d\n',rounded_IA,rounded_JA);
tissues_j)roperties_freq_rangeO;
E P S R = l;
SIGMA = 2;
%Loading Tissues' relative permittivity and conductivity vs.
%frequency values,
for frequency = 1:30
[iungs,blood,bones,heart_wali,intei'stitial_fiuid,esophagus,...
intercostal__rauscles,thoracic__wall,tumor] =
tissues_properties(frequeacy);
fprintf(l
^rintf(! .'Coordinates for Screen’®');
^rint^l
INTTX = pm i_celh + I;
E N D X = roimded_!E - INITX;
itmgs_fr(frequency,EPSR) = lungs(EPSR);
lungs_fi<fi*equency,SIGMA) = limgs(SIGMA);
biood_fr(freqx^ncy3ESR) = blood(EPSR);
blood_fr(fiequency,SIGMA) = blood(SIGMA);
bones__fr(frequency,EFSR) - bones(EPSR);
bones_fr{frequency,SiGMA) ~ bones(SIGM A);
heart_wall_:^frequency,EPSR) = heart_wall(EPSR);
heart_wall_fr(frequency,SIGMA) = hsart_wail(SIGMA);
intemtitial__fluid__fr(frequency,EFSR) = mterstitiaI_fluid(EPSR);
mterstitiaI_fiuid_fr{frequency,SIGMA) =
interstitial_fiuid(STGMA);
esophagus_fr(frequeiicy,EPSR) = esaphagus(EPSR);
esophagus_fr(frequency,SIGMA) - esophagus(SIGMA);
intercostai_muscles__fr(frequency,EPSR) ~
mtercostal_muscles(EPSR);
intercostal_muscles_fr(frequsncy,SIGMA) =
intercostal_muscles(SIGMA);
thoracic_wall_fr(frequency,EPSR) = thoracic_wal!(EPSR);
thoracic__wail_fr(frequency,SIGMA) = thoracic_waIl(SIGMA);
INITY = JC + (height__ceils/2) + dist__torso_scat_Ee!d__cells +...
dist scat field ml cells;
roimded_INITY = roimd{INITY);
rounded ENDY ==roimded__INITY;
^m ntS{ I,'INITX = %d EN DX = %d\n',INITX,ENDX);
fp rm tfil , ™ t Y = %d EN DY %d\n',roimded__INrrY/ounded_ENDY);
i^rintfi^l.Eoints to place image in Paint\n');
IC - (width_cells/2);
y f = JC + (hei^it_cells/2);
roimded_xf = round(xf);
roundedD f == round(yf);
^rintfl^ 1,’x f = %d y f = %d\n’,rounded_xCroimded_3 d);
118
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
tumor_fr(fi‘equency,EPSR) = tumor(EPSR);
tumQr__&(frequency,SIGMA) = tmnor(SIGMA);
end
thoracic_wall(COND) = skindry_cond*tw(i) -f fat_cond*tw(2) +
rauscles_cond^tw(3);
%The values o f tumor properties are obtained using a iineai- relation
o f the following tissues:
biood__prop = biood_prop_freq_rang8;
blood_epsr ~ blood_prop(FREQ,EFSR);
blood_cond ~ bloodj3rop(FREQ,COND);
bioodvessel__pfop = bloodvessel_prop_freq_range;
bloodvessel_epsr = bloodvessel_prop(FREQ,EPSR);
bloodvessel__cond - bloodvesseij3rop(FREQ,COND);
mirsc!es_prop = muscies_prop_freq_range;
muscles__epsr - muscles_j3rop(FIlEQ,EFSR);
muscles_cond = muscles_prop(FREQ,COND);
Q) “ tumorj>roportion;
tumor(EPSR) = blood_epsr^'tp(i) + bioodvessel_epsr^tp(2) +
muscies_epsr*tp(3);
tumor(COND) - Mood_cond'^p(l) bloodvessel_cond^tp(2) +
n2uscles_cond*tp(3);
%=
%This function assigns the values o f the tissues properties at a
%specific frequency given by FREQ (in GHz)
%In the 30 K 2 array tissue_prop_freq_range are contamsd the
%relative permittivity in the column 1, v/iiile the conductivity in tlie
%coiumn 2.
%The value o f FREQ is 1,2,.,.,30
function
[iimgs,biood,bones,heart_wali5mteratitial_fluid,esophagus,...
mtercostal_muscies,thoracic_wail5tumor] =
tissues_properties(FREQ);
E P S R = 1;
C O N D = 2;
% tissu e(E F S R )R ela tiv e pennittivity
% tissue(C O N r>)C onductivity S/ra.
lungs_prop = lu n ^ j5rcp_freq_range;
lm gs(EPSR )= lungs_j)rop(FREQ,EPSR);
limgs(COND) = lungsj}rop(FREQ,COND);
%This function has an array containing the
%Relati¥e pennittivity = tissue_piiop(FREQ,EPSR);
%and Conductivity = tissue_prop{FREQ,COND);
% of the tissue.
b lo o d j3rop - blood_prop_fre<L^*^g®i
bIood(EPSR) = blood_prop(FREQ,EPSR);
blood(COND) = blood_pro^FREQ,COND);
%The values o f bones properties are obtained using a linear relation
o f the following tissues:
bonecorticai_prop = bonecortical_prop_freq__range;
bonecortical_epsr = bonecorticaI_prop(FREQ,EPSK.);
bonecortical_cond = bonecortical_prop(FREQ,COND);
bonecancellousjprop = bonecancelIous_j)rop_freq_range;
bonecanceUous_epsr = bonecanceUousjpix)p(FREQ,EPSR);
bonecancellousjcond = bonecancellous__prop(FREQ,COND);
bonemarrow_prop = bonemarrow_prop_fr'eq_range;
bonemarrow_epsr = bonemarrow_prop(FREQ,EPSR);
boneraarrow_cond = bonemarrow_j5rop(FREQ,COND);
bp = bone_propoition;
bones(EPSR) = bonecortica!_epsr*bp(l) +
bonecancellous_epsr*bp(2) + ...
boneinarrow_epsr*bp(3);
bones(COND) - bonecortical__cond*b|^l)
bonecancellous_cond’*‘bp(2) + ...
bonemarrow_cond^p(3);
heart_wallj)rop = heart_wal!_j>iop_freq_raiige;
heart_wall(EPSR) = heait_wa!l_prap(FREQ,EPSR);
he^ _w ail(C O N D ) = heart_wall_prop(FREQ,COND);
interstitial_fluid_prop = int_fluid__prop_freq_range;
interstitsal_fiuid(EPSR) - interstitial_fliud_prop(FI!EQ,EPSR);
interstitial__flu!d(CONB) = iiiterstitia]__fluid_prop(FREQ,COND);
esophagusjsiop = esophagus_prop_freq_range;
esophagus(EPSR) = esophagusjprop(FREQ,EFSR);
esophagus(COND) - esopbagus_prop(FREQ,COND);
ratercostal_muscles_prop = muscies_j)rop_freq__range;
intercostal_muscles(EPSR) =
mtercostal_musc!es_prop(FREQ,EFSR);
mtercostal_muscles(COND) ="
intercoslal_m sciesj3rop(FREQ ,CO ND);
%The values o f thoracic w all properties are obtained using a linear
relation o f the following tissues:
sidndry_prop = skindry_prop_fre£|_rmige;
skindry^epsr = skindry_prop(FREQ,EPSR);
s1cmdiy__ccnd = slcmdry_prop(FREQ,COND);
fat_prop = fet_j>rop_fi^qjrange;
fat_epsr = fat_prop(FREQ,EPSR);
fet_cond = fatjprop(FREQ,COND);
muscles_prop == muscles_jjrop_freq_range;
muscies_epsr = muscles_prop(FREQ,EPSR);
muscles_cond = muscles_prop(FREQ,COND);
tw = thoracic_wallj5roportion;
thoracic_wall(EPSR) ~ skiiK!ty_epsr^tw(l) + frit_epsr*tw{2) +■
muscles_epsr*tw(3);
function [Ip] = lungs_prop_freq_range();
lp = [
21.8 2 5 ,0 .4 7 4 0 6 ; 2 0 .7 9 1 ,0 .6 8 5 2 2 2 0 .1 3 ,0 .9 6 8 9 4 19.541,1.3181;
18.966, 1.722;
18.394,2.1695; 17.8 2 3 ,2 .6 5 0 4 ; 17.256, 3.1555; 16.697, 3,6767;
16.149,4.2073;
15.614,4.7414; 15.0 9 6 ,5 .2 7 4 3 ; 14.5 9 5 ,5 .8 0 2 1 ; 1 4.113,6.322;
13.651,6.8315;
13.209,7.3288; 12.7 8 6 ,7 .8 1 2 8 ; 1 2 .384,8.2825; 12, 8.7375;
11.635, 9.1773;
1 1.288,9.602; 10.959, 10,012; 1 0 .646,10.406; 10.349,10.787;
10.066, 11.153;
9.79 8 2 ,1 1 .5 0 5 ; 9.5436, 11.844; 9.3017, 12.171; 9.0719,12.485;
8.8534, 12.787];
%'
%This function has ati array containing the
%Relative pennittivity = tissuej5rop(FREQ,EPSR);
%and Conductivity = tissue_prop(FREQ,COND);
% of the tissue.
function [bp] = blood__prop_freq_range();
bp = [
6 1 .0 6 5 ,1 .5 8 2 9 ; 5 9 .0 2 2 ,2 .1 8 6 1 ; 57.3 5 3 ,3 .0 4 9 8 ; 5 5 .677,4.1336;
5 3 .95,5.3951;
5 2 .184,6.7949; 5 0 .3 9 7 ,8 .2 9 7 2 ; 48.61, 9.8711; 46.843, 11.49;
45.109, 13.131;
43.421, 14.777; 4 1 .7 8 6 ,1 6 .4 1 2 ; 4 0 .2 1 2 ,1 8 .0 2 5 ; 38.701, 19.607;
37 .256,21.152;
35.878, 22.654; 34.566, 24.111; 3 3 .3 1 8 ,2 5 .5 1 9 ; 3 2 .134,25.879;
31.011,28.189;
29.9 4 6 ,2 9 .4 5 1 ; 2 8 .9 3 7 ,3 0 .6 6 4 ; 2 7 .9 8 2 ,3 1 .8 3 1 : 27.0 7 7 ,3 2 .9 5 1 ;
26.219, 34.027;
25.4 0 7 ,3 5 .0 6 1 ; 24.637, 36.053; 23.907, 37.006; 23.215, 37.922;
22 .559,38,801];
%This function has an array containing tiie
%Relative permittivity ~ tissuejjrop(FREQ,EPSR);
%and Conductivity = tissue_prop(FREQ,COND);
% of the tissue.
function [bcp] = honecorticaI_prop_freq_rangeO;
119
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
%From: The Biomedical Engineering Handbook, 2ed Vol. I,
%CRC Press, IEEE Press, 2000.
% bp{2).- Corresponding to bonecancelious
%bp(3) .- Corresponding to bonemarrow
%Soft material filling a central cavity, especially o f bones.
%Soft tissue in die cavities o f bones.
%From: The Harper Collins Illustrated M edical Dictionary
%Harper Resource, HArper Collins Publishers, 4ed.. 2001
%The proportion w as obtained ftom figures from the above
%rr:-entioned books.
%And fromhttp;//www.bartleby.com /!07/28.htm l
%Structure. The ribs consist o f highly vascular cancellous
%tissue, enclosed in a thin layer o f compact bone.
bsp = [
12.363.0.15566; 1L654, 0.31007; 11.066,0.50624; 10.532,
0.72752; 10.04, 0.96228;
9 .5869,1.2026; 9 .171,1.4431; 8.7896, 1.6801; 8.4401, 1.9115;
8.1197, 2.1359;
7 .8 2 6 ,2 .3 5 2 3 ; 7 .5563,2.5605; 7 .3 084,2.7604; 7.0802,2.9519;
6.8598,3.1355;
6.6755, 3.3113; 6.4957, 3.4797; 6 .3 29,3.6411; 6.1741,3 .7 9 5 8 ;
6.0301, 3.9441;
5.8959,4.0865; 5.7706,4.2233; 5.6534,4.3547; 5.5436,4 .4 8 1 2 ;
5.4406, 4.603;
5 .3 439,4.7203 ; 5.2529,4.8335; 5.1671,4.9427; 5.0862,5.0482;
5.0097,5.1502];
b p ( l ) - 0 .3 ; bp(2) = 0.6; b p { 3 ) - 0 .1 ;
%This function
an array containing the
%Relative permittivity = tissuejprop{FREQ,EPSR);
%and Conductivity ~ tissue j>rop(FREQ,COND);
% of the tissue.
%This function has an array contaming the
%Relative permittivity = tissue_prop(FREQ,EPSR);
%and Conductivity = tissue_j)rop{FREQ,COND);
% of the tissue.
function [bcp] = bonecanceiious J3rop_^freq_range0;
function [hwp] = heart_wall_pmp_freqj:angeO;
bcp = [
20.584.0.36395; 19.086,0.6522; 17.943,1.0062; 16.946,1.3988;
16.05, 1.8116;
15.238,2.2318; 14.502,2.6512; 13.831, 3.0643; 13.22, 3.4676;
12.661,3.8591;
12.15,4.2376; 11.682, 4.6024; 11.252, 4.9536; 10.856,5.2912;
10.491, 5.6157;
10,154, 5.9275; 9.841, 6.2271; 9 .5509,6.5151; 9.2811,6.7 9 2 1 ;
9.0297, 7.0586;
8.7951,7.3152; 8.5757,7.5625; 8.3702,7.8009; 8.1773,8.0309;
7.9962,8.2531;
7.8256, 8.4677; 7.6649, 8.6754; 7.5132, 8.8763; 7.3698, 9.0709;
7.2341,9.2595];
hwp = [
59.2 9 ,1 .2 8 3 6 ; 5 5 .8 1 6 ,1 .9 1 1 7 ;
50.274, 4.8626;
4 8 .6 2 ,6 .1 2 2 3 ; 4 6 .9 8 6 ,7 .4 7 2 9 ;
42.241, 11.835;
40.735, 13.329; 39.277, 14.819;
3 5 .222,19.171;
33.9 8 2 ,2 0 .5 6 1 ; 32.799, 21.913;
29.576, 25.724;
2 8.6 0 6 ,2 6 .9 1 ; 2 7 .6 8 4 ,2 8 .0 5 5 ;
25.188,31.242;
24.4 3 9 ,3 2 .2 2 7 ; 2 3 .7 2 7 ,3 3 .1 7 4 ;
2 1 .798,35.808];
%This function has an array containing the
VoRelative permittivity - tissue_prop(FREQ,EFSR);
“/oand Conductivity = tissue_prop(FREQ,COND);
%of the tissue.
5 3 .7 3 7 ,2 .7 2 8 8 ; 5 1 .961,3.7216;
4 5 .3 7 4 ,8 .8 8 9 4 ; 43.79, 10.35;
37.87, 16.295; 36.518, 17.747;
3 1 .6 7 1 ,2 3 .2 2 5 ; 3 0 .597,24.496;
2 6 .8 0 9 ,2 9 .1 5 7 ; 2 5 .9 7 8 ,3 0 .2 2 ;
2 3 .0 5 1 ,3 4 .0 8 6 :2 2 .4 0 9 ,3 4 .9 6 3 ;
%This fimction has an array contaming the
%Relative pennittivity = tissue_prop(FREQ,EPSR);
%and Conductivity = tissue_prop(FREQ,COND);
% of the tissue.
function [bmp] - bone!mrrow_prop_freq_rangeO;
fimction [ifp] = int_fluid_prop_freq_range();
bmp = [
5.4854.0.042803; 5.3478,0.07615; 5.2378,0.12085; 5.1356,
0.17413; 5.0379, 0.23379;
4.9441.0.29808; 4.8541,0.36563; 4.7679, 0.43532; 4.6856,
0.5063; 4.607,0.57788;
4.5322, 0.64952; 4.461, 0.7208; 4.3933, 0.7914; 4.329, 0.86107;
4.268, 0.92961;
4.2 1.0.9969; 4.1549, 1.0528; 4.1025, 1.1273; 4.0528, 1.1904;
4.0055, 1.252;
3.9606,1.312;
3,9178, 1.3706; 3.8771, 1.4278; 3.8383, 1.4834;
3.8014, 1.5377;
3.7661, 1.5906; 3.7325,1.6421; 3.7003, 1.6923; 3.6697, 1.7413;
3.6403, 1.789];
6 8 .8 7 5 ,1 .6 6 7 3 ; 6 8 .4 7 2 ,2 .1 5 5 6 ; 6 7 .8 1 7 ,2 .9 5 5 9 ; 6 6 .9 2 3 ,4 .0 4 9 4 ;
65.8 1 ,5 .4 1 1 ;
64.5 0 2 ,7 .0 1 1 2 ; 6 3 .0 2 6 ,8 .8 1 7 2 ; 6 1 .4 1 1 ,1 0 .7 9 5 : 59.683, 12.909;
57.872, 15.126;
56.002, 17.414; 54.098, 19.745; 5 2 .1 8 ,2 2 .0 9 3 ; 50.2 6 7 ,2 4 .4 3 4 ;
48.374, 26.75;
46.516, 29.025; 44.701, 31.246; 4 2 .9 3 8 ,3 3 .4 0 4 ; 4 1 .2 3 3 ,3 5 .4 9 1 ;
39.59, 37.501;
3 8 .013,39.432; 3 6 .5 0 2 ,4 1 .2 8 1 ; 35.0 5 8 ,4 3 .0 4 8 ; 33.6 8 1 ,4 4 .7 3 3 ;
32.37, 46.338;
31.123, 47.865; 2 9 .9 3 8 ,4 9 .3 1 5 ; 2 8 .8 1 3 ,5 0 .6 9 2 ; 27.7 4 5 ,5 1 .9 9 9 ;
26.733, 53.238];
%This hinction assign the proportion, in which the values o f the
%properties o f the tissues;
% bonecortical, bonecancellous, and bonerrmrrow
% are used to build "bone".
% The sum o f them must be i.O
%This fimction has an array containiiig the
%Relative permittivity = tissuejprop(FREQ,EPSR);
%and Conductivity = tissue_prop(FREQ,COND);
% of the tissue.
function [bp] = bonej3ropoition()
fimction [ep] = esophagus_prop_freq_rangeO;
% bp(l) Corresponding to bonecortical
%The dense compact bone found throughout the shafts o f long bones
%such as the femur, tibia, etc. A lso found in the outer portions
% of other bones in the body.
ep = [
6 4 .797,1.2316; 6 2 .8 9 2 ,1 .8 4 3 5 ; 61.2 6 8 ,2 .7 2 9 2 ; 59.6 1 1 ,3 .8 4 7 6 ;
57.89, 5.1565;
120
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
36.118,6.6165; S 4 .3 i5 , 8.1918;
48.923, 13.314;
47.184, 15.075; 45.491, 16,834;
40,757,21.978;
3 9 .304,23.623 ; 37.914, 25.223;
34.123, 29.735;
3 2 .979,31.14; 3 1 .892,32.495;
28.946, 36.269;
2 8 .061,37.434 ; 27.2 2 ,3 8 .5 5 5 ;
24.94,41.672];
52.503,9.8508; 50.7,11.566;
3.9581,
3,7993,
3.7641,
3.6385,
43.853, 18.577; 42.274,20.295;
1.3181; 3.9155, 1.3766; 3.8748, 1.4336; 3.8361, 1.4893;
1.5434;
1.5962; 3.7305, 1.6477; 3.6985, 1.6979; 3.6678, 1.7467;
1.7944];
36.588,26.776; 35.325,28.28;
30.859,33.8;
29.878,35.058;
%Tiiis function has an array contaming the
% Rsiativs permittivitj' = tissue_prop(i7REQ,EPSR);
%and Conductivity = tissue_prop(FREQ,COND);
% of the tissue.
26.421,39.634; 25.662,40.673;
fimction [it^)] = museles_prap_freq_rangeO;
%TUs fimction has an array containing the
%Reiative permittivify = tisstte_prop{FEEQ,EPSR);
%and Conductivity = tissue_prop{FREQ,COND);
% of the tissue.
mp = [
54.8 1 1 ,0 ,9 7 8 1 9 ; 53.29, 1.4538; 52.0 5 8 ,2 .1 4 2 1 ; 50.321,3.0155;
4 9 .5 4 ,4 .0 4 4 8 ; 4 8 .2 1 7 ,5 .2 0 1 9 ; 4 6 .8 6 5 ,6 .4 6 0 7 ; 4 5 .497,7.7978;
4 4 .1 2 6 ,9 .1 9 2 3 ; 4 2 .7 6 4 ,1 0 .6 2 6 ; 41.419, 12.083 ; 40.1 0 1 ,1 3 .5 5 ;
38 .8 1 4 ,1 5 .0 1 6 ; 37.564, 16.472; 36.354, 17.911; 35.186, 19.326;
34 .0 6 1 ,2 0 .7 1 2 ; 3 2 .9 8 ,2 2 .0 6 7 ; 3 1 .9 4 4 ,2 3 .3 8 8 ; 3 0 .951,24.673;
30 ,2 5 .9 2 1 ;
2 9 .0 9 2 ,2 7 .1 3 1 ; 2 8 .2 2 4 ,2 8 .3 0 3 ; 27.395,29.437;
26 .6 0 3 ,3 0 .5 3 4 ; 2 5 .8 4 8 ,3 1 .5 9 5 ; 2 5 .1 2 7 ,3 2 .6 1 9 ; 24.4 4 ,3 3 .6 0 9 ;
23 .7 8 3 ,3 4 .5 6 5 ; 2 3 .1 5 7 ,3 5 .4 8 7 ];
function [mp] = muscles_prop_freq_rangeO;
mp = [
5 4.811.0.97819; 53.29, 1.4538; 52.058,2.1421; 50.821, 3.0155;
49.54, 4.0448;
48.217,5.2019; 46.8 6 5 ,6 .4 6 0 7 ; 45.497,7.7978; 44.126,9.1923;
42.764, 10.626;
41.419, 12.083 ; 40.101, 13.55; 38.814, 15.016; 37.564, 16.472;
36.354, 17.911;
35.186, 19.326; 3 4 .061,20.712; 32.98,22.067; 31.944,23.388;
30.951,24.673;
3 0 ,2 5.921;
2 9 .092,27.131; 28.224,28.303; 27.395,29.437;
26.603, 30.534;
25,848,31.595; 2 5 .127,32,619; 24.44, 33.609; 23.783, 34.565;
23.157, 35.487];
%
“/©This function assign the proportion, in which the values o f the
%propeities o f the tissues:
% s k ^ r y , fat, and muscle
% are used to build "thoracic wall".
% The sum o f them must be 1.0
fimction [tw] = thoracic_waIijproportionO
% tw (l)C o r r e sp o n d in g to skindry
% tw (2 )C o rresp o n d in g to fat
% tw (3 )C o rresp o n d in g to muscle
%This function has an array containing the
%Relative permittivity = tissue_prop(FREQ,EPSR);
%and Conductivity = tissue_prop(FREQ,COND);
% of the tissue.
%The proportion w as obtained by checking the image 1451 from
%the Thorax from the Visible Human Project.
t w ( ! ) - O .I ; tw(2) = 0.5; tw (3) = 0.4;
function [sdp] = skindryjrop_freq_rangeO;
sdp = [
40.936.0.89977; 38.568,1.2654; 37.45, 1.7406;
35.774, 3.0608;
34.946,3.8912; 34.084,4.8175; 33.184,5.8242;
31.29, 8.0138;
30.313,9.1658; 29.327, 10.337; 28.342, 11.514;
26.401, 13.847;
25.458, 14.986; 24.539, 16.097; 23.649, 17.176;
21.962,19.224;
2 1 .168,20.19; 2 0 .409,21.114; 19.684,21.998;
18.335, 23.644;
17.709,24.408; 17,115,25.134; 16.552,25.824;
15.51,27.099];
36.587,2.3404;
%This fimction has an array containing the
%Relative permittivity - tissue_prop(FREQ,EFSR);
%and Conductivity = tissuejprop(FREQ,COND);
% of the tissue.
32.25,6.8949;
27.364, 12.688;
function [bvp] = bloodvessel_prop_freq_™>tS®0;
22.789, 18.219;
bvp = [
44 .5 6 1 ,0 .7 2 8 6 6 ; 4 3 .0 8 9 ,1 .1 7 0 8 ; 4 1 .8 5 5 ,1 .8 0 8 ; 40.5 9 7 ,2 .6 0 6 6 ;
39.2 9 5 ,3 .5 3 3 ;
3 7 .9 6 3 ,4 .5 5 6 2 ; 3 6 .6 1 9 .5 .6 4 8 9 ; 3 5 .2 8 ,6 .7 8 7 4 ; 3 3 .961,7.9522;
32.6 7 3 ,9 .1 2 6 8 ;
31.425, 10.298; 30.223, 11.456; 29.07, 12.593; 27.97, 13.703;
26.922, 14.782;
25.926, 15.827; 24.982, 16.835; 24.089, 17.807; 23.244, 18.742;
22.445, 19.64;
2 1 .6 9 1 ,2 0 .5 0 1 ; 2 0 .9 7 8 ,2 1 .3 2 8 ; 2 0 .3 0 5 ,2 2 .1 2 ; 19.67,22.879;
19.069,23.606;
1 8 .502,24.303; 17.966,24.971; 17.458,25.611; 16.978,26.225;
16.524,26.814];
18.993,22.841;
16.017,26.478;
%This function has an array containing the
%Relative pennittiviSy - tissue_prop(FREQ,E^S^);
%and Conductivity = tissue_j>rop(FREQ,COND);
% of the tissue.
function [fp] = fat_prop_freq_range();
K
ip=[
5.447.0.053502; 5.3276,0.085918; 5.2239,0.13004; 5.1249,
0.1829; 5.0291.0.24222;
4.9367.0.30623; 4.8476,0.37353; 4.7622,0.44301; 4.6804,
0.5138; 4.6023, 0.58521;
4.5278, 0.65669; 4.457, 0.72783; 4.3896, 0.79829; 4.3255,
0.86783; 4.2647, 0.93625;
4.2068, 1.0034; 4.1519, 1.0693; 4.0997, 1.1337; 4.0501, 1,1966;
4.003, 1.2581;
%This fimction assign the proportion, in which the values o f the
%properties o f the tissues:
% blood, blood vessel, and naisd e
% are used to build "turoor".
% The sum o f (hem must be 1.0
fimction [^J = tumor_j)roportionO;
121
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
% tp (l)C o r r e sp o n d in g to blood
% tp (2 )C o iresp o n d in g to blood vessel
% tp (3 )C o r resp o n d in g to muscle
else
i^rint^fp,'movie:
end
%The proportion w as obtained according to the values used in
%diifereEt papers, mainly;
%"Microv/ave Imaging for Breast Tumor: 2D Forward and Inverse
Methods."
:^rinti(^,*movie_period:
%d\n’5movie_period);
tjrintf(ip,'ITERATIONS.-Vi’);
%rintf(^,'iisteps:
%d\n’,nsteps);
%
fjprmtf(^,’ax!s:
^rintl(:^/initx:
no‘\n');
%dW,axis);
%d^',INITX);
%d\n’,INITY);
^irintfi^^^/endx:
%d\n'jENDX);
ferintf(^,'endy:
% dW ,ENDy);
inntfi;tj;PROPERTIES„FILES.-\n');
^rintf(4>,'GA: Ga\n');
tp ([) = 0.60;tp(2) “ 0.20; tp{3) = 0.20;
%Generate__dimensions_hO
%This fimction generate the file dimensions.h
GbXn*);
function Generate_diir!ensions_h{IE,JE,iC,JC,IA,JA)
iprintf(Q);STEADY_STATS.-\B');
^nntf(^,'Last_j)eriod: 250W);
t.rmt^^.'PROPERTIES_FILE_I0.-'\n’);
^ rm tf(^ ,’Prop__ID: prop_ID\n');
$rintfi;^;AVERAGE_STEADY_STATE.>\ii');
^rin.tf(:^,’slart_averaging_step: 1000\n’);
fclose(^ );
^ = fopen('dimensions.h’,'wt’);
^rintf(:^/dimensioi2S.li\n');
^rintfi^Q3,'\n');
§)rintf(^,V* Dimensions ofF D T D space *An');
^rint^%,'#define IE %d\n’,IE);
fprintfl[fp,'#define JE %dV,JE);
43rintf(^,V’®Total/scattered field boundaries *An');
^srintfi^fp.'Mefme lA ^d\n',IA);
^ rintf(^,'#defm e IB (IE-IA-1 )\n’);
^rintfl:Q>.'#define JA %dW,JA);
Q)rintf(^,’#define JB (JE-JA-l)\n’);
Coordinates ofF D T D space center ®An');
^rintf(fp,*#define IC IE/2W);
Q5rintf(f^;#defme JC JE/2\n');
fciose(f|s);
%'
%Function Generate_input_txt()
%This fimction generate the file input.txt
function Generate_input_txt(deita_x,source_frequency,tO,spread,..
npml,movie 5movie_period,nsteps,...
axis,IN ITX ,IN ITY ,EN D X ,EN D ^
^ = fopen('mput,txt','wt’);
This file contains the parameters for the program\n');
2dfdtd.c.\ri');
^rintfi^fp,’* B y the time, the file has a rigid format, that nieans,\n’);
^rint^fp,'* the format must not be changed.\n');
:^rintfl^fp,'* It is highly recommended to make a backup o f this\n');
file, and then work with the copy.\n');
^rintf(i|3,'* There must be a least o i^ space between the name o f
the\n');
^rintfi^fp,’* variable and its value.'m');
^rm tf(^,'^ Writing in this file is allowed, but only above the\n');
following line.\n');
^rintf(^;START_PAIlAM ETERSW );
4>rintfi;^,DIM ENSI0NS.-V);
^rintf(;^,'deita_x;
%l0.9fvn',delta_x);
^rintf(fp,*SOURCE.-\n');
fprintf(fi3,'source__frequency: %d.e9\n’,source_frequency);
^rintf(i,'tO:
%4.2f\n',t0);
^rintfi'^,'spread;
%4.2f\n’,spread):
^rintfi:ip,'PML,*\n');
%d\n’,npml);
^rintf(§),'MOVIE.-\n’);
i f (movie— !)
^rintf(^,'movie:
yesW);
122
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX D
MATLAB PROGRAM “IMAGE2CONSTANTS”
123
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
%p2operties are according to the frequency o f interest.
%Program: miage2constants
%Using as input the inmge o f the biologicaJ phantom, the anuys
%contaimng the constants Ga and Gb are generated.
clear all;
cic;
switch (freq)
case 1
[phantom,map] = im readCtotso_0iG Hz_I58xl58.bm p',bmp');
[pha,ntoni,map] =
imread('torso_01G H z_l 58x 15S_p2_d3 .bn^^bmp');
- fopen('prop_ID_l_l58_p2_d3'j*wt');
= fopen(’tissuesjproperties_iG H z.h’,'wt');
^ 3 = fopen(’Permittivities_l_!58_p2__d3',’vk3;
case 2
[phantom,map] = imreadCtorso_02GH2_284x284.bmp',’bmp');
[phantom,map] =
ixnread('torso_02GHz__284x284_p2_d3.bmp',‘bmp');
~ fopen('prop_ID_2_2S4_p2_d3Vwt“);
= fopenC'tissuesjjroperties^lGHz.h'j’wt*);
fp3 = fopen('Permittivities__2_284__p2_d3','wt3;
case 3
[phantom,niap] = imread(’torso_03GHz_4!0x4i0.biBpVbmp');
[phantom,map] =
imread('torso__03GHz_410x410_p2_d3.bmp',brap');
fr)l = fopenCprop__ID_3_410_p2_d3Vwt’);
= fopen('tissues_j)roperties_3GHz.h','wt’);
= fopen{'Permittivities_3__410_p2_d3’,'wt’);
fonnatlong;
%Vaiues o f the different colors
[c01,c02,c03,c04,cQ5,c06,c07,cOS,c09,ciO] = colorsOitoIQ;
[c 1 i ,c 12,c 13 ,c 14,cl 5,c 16 ,c! 7,c 18,c 19,c20j - colors 11 to20;
%Associate a tissue to a color
lfeespace__col = cOl;
limgs__col = c02;
blood__coi - c03;
boBes_coi = c04;
heait_wall__col = c05;
interstitial__fluid_col = c06;
esophagus_col = c07;
intercostal_muscies_col = c08;
thoracic_wal!_coi = clO;
tumor_co1 = cl9 ;
[phantoin,map] = imread('toi^o_04GHz_532x532.brrp',bnp');
[phantom,map] =
irriread(’torso_04GHz_532x532_j)2_d3.bmp',’bmp’);
= fopenCprop_ID_4_532_p2_d3’,*wt');
= fopen(’tissues_properties_4GHz.h',*wt’);
fr)3 = fopen('Perrrattivities_4_532_p2_d3','wt’);
case 5
[phantom,map] = inimad(’torso_05GHz_652x652.bmp*,brnp3;
[phantom,map] =
imread(torso_05GHz_652x652jp2__d3.bfiip',bmp');
= fopei^’prop_ID_5_652_p2_d3Vwt');
= fopen('tissuesj5roperties_5GHz.hVwt');
§)3 = fopen(’Permittivities__5_652_p2_d3','v^');
case 6
[phantom,imp] = iraread('torso_06GHz_770x770.bnq>',biip*);
[phantom,map] =
!mrea£^’torso_06GHz_770x770Jp2__d3 .bmp'sbmp');
§)1 = fopenCprop_ID_6_770_p2_d3','wt’);
= fopen(tissues_j>roperties_6GHz.h’,'wt');
^ 3 - fopen(*Permittivities_6_770_p2_d3Vwt');
case 7
[phantom ,m ^] = miicad('torso_7GHz_i 114x1114.bmp',bmp');
^ 1 = fopenCprop_ID_7_ 1114*,Vt');
^ 3 = fopen('Permittivities_7_l 1 i4Vwt');
case 8
[phantom,map] = imread(torso_8GHz_1252xl252.bn:^Vbmp');
^ 1 = fopenCprop_ID_8_I252',’wt*);
:^3 - fopen('Permi£tivities_8_i252',Svt');
case 9
[phantom,map] = imread(‘torso_9GHz__l384xl384.bmp',bmp');
§ ji = fopen(‘prop_ID_9_1384','wt');
^ 3 -fopen('Pennittivides_9_1384',S¥t');
case 10
[phantom,niap3 =
imread('torso_IOGHz_l 196x1196_j32_d3.bmp',bmp');
fp l = fopen('propJD_IO_! !96_p 2 _ d 3 ;V f);
= fop€n(tissues_j)roper£ies_10GHz.hVwt');
^ 3 = fopen(Term ittivities_10_l !9 6 _ p l_ d r,V t');
case i 1
[phantom,map] =
imread(’to rso _ !l G H z_l 631 x l 6 3 1.bmp'jbmp');
=foper6Cprop_ID_II_i632‘,Vt');
%3 = fopen('Permittivities_il__1632Vwt');
case 12
[phantom,map] =
imread('torso_12GHz_l746x1745.bmp’,bm p’);
= fopenCprop_ID_12_i746Vwt^;
fr}3 = fopen('Permittivities_12_1746','wt');
case 13
%Tissues' IDs.
NUM BER_OF_TISSUES = 10;
FREE_SPACE = 0;
LUNGS = 1;
BLOOD = 2;
BONES = 3;
HEART_WALL = 4;
INT_FLUID = 5;
ESOPHAGUS = 6;
INT_MUSCLES - 7;
TH0RACIC_W ALL = 8;
TUM OR = 9;
%EM constants
^ ;|c 4c $ ;t> ; j c 1$ ;« ^ :jt 9|c:)(ili ;Jc tj! 4 t 4 : « # ^ i tc ^ :js 4 i ^ ^ 9^ ^ :j> $ 4 : $ 4s ^ 1^
EFSZ = 8.85419e-12;
^
ij( ^ :i: ^
C -3 e8 ;
%FDTD parameters
%*
%whiie r e s p = 0
% ddx = ii^utCEnter ddx from torso_dimensions_in_celis: ’);
% Q>rtntfl;i,'Is correct ddx - % 10.9f^ ',ddx);
% resp = input('l -> yes, 0 -> n o :');
%end
%This value o f ddx corresponds to the cell size according to
%torso_dimensions_m__cells
%It may be that later I would use ddx corresponding to 30 GHz
%or less, although the incident si^ ial could less. With this w e could
%analyze the response to a certain frequency but with a smoother
%structure.
%dt = d d x/(2*C );
%Load information o f permittivity and conductivity o f tissues for
%a specific frequency.
% tis s u e (l)= relative permittivity
% tissue{2)= conductivity
freespace(l) = 1.0;
freespace(2) = 0.0;
freq = input(Trequency to analyze? l-> 3 0 GHz: ’);
[iui^,blood,bones,heait_wali,interstitiai_fluid,esophagus,..,
intercostal_muscles,tboracic_wail,tumor] = tissues_properties(freq);
%Read image
%These files are the size according to the frequency, so am the
%propetties.
%Should w e want to use a "high-frequency*' size structure for
%a "low-frequency" analysis, w e have to m k e sure the
124
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
case 27
[phantom,map] =
irnrea<!(torsQ__27GHz_2S96x2S96.bmp',’bmp');
=ftipen('prop__iD__27___2896’,W ) ;
^ 3 = fopen('Pemiittivities_27_2896Vwt');
case 28
[phantoiiLmap] =
imread('torso_28GHz_2336x2336.buip','bmp');
[phantom,map] =
miread('torso 28GH z_2336x2336.bm p’,‘bn^')^
fopen{'prq>_ID_28_2336_pO_dOV'wtO;
= fopen('tissuesj5roperties__28GHz.h',’wt');
^ 3 = fopen('Fennittivities_2S_2336Vwt');
case 29
[phantom,rmp] =
!mread(’torso__29GHz_2994x2994.bmp’,'bi'!ip’);
= fopen{’prop__ID_29_2994',’wt');
^ 3 = fopen(TPermittivities__29_2994',Vt*);
case 30
[phaBtom,map] =
imread{’torso_30GHz_3040x3040.bnipVbmp');
Q?! = fopen(’prop_ID_30__3040','T^');
= fopen(Fermittivities_30_3040Vwt‘);
otherwise
disp(W rong frequency!*);
end
[pliantom,roap] =
imread('torso__13GHz_l 856x1 S56.bmp‘,’btnp');
~ fopen(’prop_ID_l3__iS56’,'wt’);
^ 3 ~ fopen('PermitUvities_13__i856','wf);
case 14
[pbantoni,mapl =
imreadCtorso_14GHz_l 552x1 S52.bmpVbn^');
[phantom,imp] =
imread(’torso_I4G H 2_1552s!552_p2_d3.bm p',’b5i^');
Ipl = fopen('propJD_14_i552_p2_d3Vwt');
“ fopen(’tissues_properties_14GH2.hVwt');
^ 3 = fopenCPermittivities_l4_l552','wt');
case 15
[phantom,map] =
ioiread('torso__15GHz_205 8x2058.bmp*,brrip’);
“ fopen('prop__ID__l5__2058’,'^ ');
f^3 = fopeE(‘Permittivities_i5_205SVwt');
case 16
[phantom,niap] =
mmBad(’torso_l 5GH z_2150x2 i 50.bmp','bmp’);
=fopen('prop_ID _!6_2l50V w t’);
^ 3 ~ fopen('Perraittivities__16_2!50','wt');
c ^ e 17
[phaiitom,map] =
iimead(torso__l 7GH2__223 8x223 8,bmp’,l3rnp’);
Ipl = fopenCprap_ID_l7_2238',‘wt’);
^ 3 = fopenCPennittivities_17__223S',‘wt');
case 18
[phantom,map] =
iraread('torso_l 8G H z_i S40x! 840.bmp',’b itp ’);
[phantom ,m^] =
imread(’torso_l 8G H z_l 840x 1840_p2_d3 .brapVbmp');
= fopenCprop_ID_18__1840_j)2_d3','wt');
= fopen('tissuesjpropeities_i8GHz.hVwt');
fp3 - fopen('Pennittivities_!8_1840','wt');
case 19
[phantoni,map] =
imread('torso_i9GHz_240Gx2400.bn^’,’bmp');
^ 1 =fopenCprop__ID_19_2400Vwf);
^ 3 = fopen('Permittivities_19_2400','wt');
case 20
[phantom,map] =
imread(torso_20GHz__2474x2474.bmp','bn::p');
^ 1 = fop0n('prop__ID_2O_2474Vwt*);
^53 = fopen('Perraittivities__20__2474Vwt');
case 21
[phantom,map3 =
iETiread('torso_21GHz_2544x2544,bmp','binp’);
§)1 = fopen('prop_ID _21_2544'>t’);
^ 3 = fopen('Permittivities_21_2544Vwt');
case 22
[pfaantom,map] =
imread('torso_22GHz__2610x2610.bn^',‘bmp');
Ijsl - fopei^’prop_ID_22_26iO’,‘w t‘);
1^3 = fopen('Permittivities_22_2610',’wt');
case 23
[phantom,mapJ =
!mread('torso__23GHz_2674x2674.bmp','brr43');
=fopen('prop__lD__23_2674',\vt');
%3 = fopen('ipermittivities_23_2674','wE');
case 24
[phantom,map3 =
imreadCtorso_24GHz_2734x2734.bn^','bmp');
ip l = fopen('prop_ID_24__2734Vwt’);
^ 3 “ fopen('Penmttivities_24_2734'5'wt’);
case 25
[phantom,raap] =
3mread('torso_25GHz_2790x2790.bnip’,bn:^');
1^1 = fopenCprop_ID_25_2790’,V f);
fp3 = fopen(’Fermittivities__25_2790',*wt');
case 26
[phantom,map] ==
imread('torso_26GHz_2844x2844.bmp',l3mp');
^ 1 - fopen('prop_ID_26_2844Vwt*);
^ 3 = fopen('Permittivities_26_2844’j‘wt');
%Start time
time! = f!x(clock);
^5rintf(i,'Initial time; %d:%d:%d\n',timel(4),time!(5),tiraei(6));
[lEGE,layers]=size(phantom);
% A ssi^ m g ones to the array mask,
mask = oiks(IE,JE);
%Scanning the image
fo r j= l JE
fori= l:IE
tissue = phantom(iJ,:);
if (tissue — freespace_col)
epsr = freespace(i);
sigma = freespace(2);
ID = FREE_SPACE;
Emsk(i,i)“ 0;
else if (tissue ~ lungs_col)
epsr = hmgs(l);
ID = LUNGS;
else if (tissue = = blood_col)
epsr = blood(l);
ID -B L O O D ;
else if (tissue = bones_coI)
epsr = bones(l);
ID = BONES;
else if (tissue — heait_wall_col)
epsr = heart w all(l);
ED = HEAilT_WALL;
elseif (tissue ~ mt€TStitial__fluid_col)
ep sr= mteratitiai_fluid(l);
ID = INT_FLUID;
else if (tissue = esophagus_co!)
epsr = esophagus(l);
ID = ESOPHAGUS;
e b e if (tissue = intercosta!_muscies_col)
ep sr= intercostai_muscies(i);
ID - INT_MUSCLES;
else if (tissue = thoracic_waii_col)
epsr = thoracic_wali(l);
ID - THORACIC_WALL;
else if (tissue = = tumor_col)
epsr = tumor(l);
ID = TUMOR;
else
125
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
%7.5fn',epsr);
fprktfC^t.Wdefme BLOOD_EPSR
t>rint^it,'#deflne BLOO D_SiG M A
%7.5Te',sigma);
epsr = bones(l);
s ig n a = bones(2);
%7.5fn',8psr);
§)rmtf(%t,'#defme BONES_EPSR
%7.5f\n',sigma);
^ r m tf ftit.Wdefine BONES_SIGM A
epsr = heart_wall{l);
sigma = heart_wall(2);
§5rmtfl:#»t,'#define HEART_W AIX_EPSR
%7.5f«i',epsr);
^rintf(43t,'#defme HEART_W ALL_SIGMA
%7.5f(n',sigma);
spsr = interstitial_fluid(l);
sigma = interstitial_fluid(2);
%7.5f\n’,epsr);
^)rmti(ipt,’#define INT_FLUID_EPSR
%7.5ftn',sigma);
t5rintf({i)t,'#defme INT_FLUID_SIGMA
epsr = esophagus(l);
sigma = esophagus(2);
%7.5f\n',epsr);
if>rintf{%t,'#defme ESOPHAGUS__EPSR
§>rmti)^t,'#define ESOPHAGUS_SIGMA
epsr - inten;ostal_muscles(!);
sigma = mtercostal_muscles(2);
%7.5f\n’,epsr);
§5rintl(fi)t,'#define 1NT__MUSCLES_EPSR
iprijitfl'it,’#define INT_M USCLES_SIGM A % 7.5f<sigm a);
epsr = thoracic_wall(l);
sigma = thoracic_wall(2);
t>rintfl:ft)t,'#defma THORACIC_W ALL_EPSR %7.5f\n',epsr);
^rintf(4)t,'#defme THO RACIC_W AIX_SIGM A %7.5f\ti',sigma);
epsr = tumor(l);
sigma = tumor(2);
^rintfl:^t,'#defme TUMOR_EPSR
%7.5fti',epsr);
irm tf(4)t.'#define TUM OR_SIGM A
%7.5frf,sigma);
iprintf(l,XJiikiiown Tissue!%’);
:^rintl(l /A t position %2,0f,%2.0f\n',ij);
Taint cocrdmales; %2.0f,%2.0fin'J--l,i-l);
dummy = input('Yeriiy the in^gs');
end
IS5rinta:ipl,'%d MD);
^ r in tf (i3 /% i0 .9 f’,epsr);
end
lprintf(§5l,‘\n');
end
fclose(:^!); fclose(^3);
%7.5f\ii'5sigma);
%Saving the mask to use it later in some postprocessing.
%save mask_torso mask;
%clear mask;
%Show image
f!gure(i);
imagesc(phantom);
colormap(imp);
axis square;
xIabelCx');
ylabel(y>;
clear pliantom;
load mask_torso;
figure(2);
surf(mask);
xlabelCx');
ylabel('y');
axis([l,IE,l,JE]);
axis square;
view(90,90);
shading flat;
colormap(piiik);
Y E S -1 ;
PROP ID = NO;
fciose(lpt);
%Start time
^rintfl^lslnitial time: %d:%d:%dW,tmiel{4),timei(5),tin:Kl(6));
%End time
time2 = fix(clock);
fyrim^l,'Etid time: %d:%d:%d\n',time2(4),time2(5),time2(6));
N O -0 ;
% A ssign numerical values for the first 10 colors.
iflTRO P.ID)
strlD - Cprop_ID__10_1196_p2_d2');
load(strlD);
figure(20);
surf(prop_ID_l 0 _ 1 196_p2_d 1);
xlabelCx');
ylabelCy*);
axis([U E ,l,JE ]);
axis square;
view(90,90);
shading flat;
coiorbar;
end
fimction [c01,c02,c03,c04,c05,c06,c07,c08,c09,cl0] =
colorsOltolOO
c 0 1 (i,1 ,:) = [255,255,255];
c 0 2 ( l,l,; ) - [192,192,192];
c 0 3 (l,i,:) = [255,0,0];
c 0 4 (!,l,:) = [255,255,0];
c 0 5 ( l,l,: ) = [0,255,0];
c 0 6 (l,l,:) = [0,255,255];
c 0 7 ( l,l,; ) = [0,0,255];
c 0 8 ( l,l,: ) = [255,0,255];
c 0 9 ( l,l,: ) = [255,255,128];
c ! 0 ( l ,l,: ) = [0,255,128];
%Generating tissues_properties.h
fimction [ c U ,c l2 ,c ! 3 ,c l4 ,c l5 ,c ! 6 ,c ! 7 ,o l8 ,c l9 ,c 2 0 ] =
colors! lto20()
c l2 (i,l,
[128,128,255];
c l 1(1,1,:) = [128,255,255];
cl4(1,l,
[255.128.64];
[255,0,128];
c l 3 ( l ,l ,
[64,0,255];
[128,64,0];
c l 5 ( l ,l ,
c l6 ( l,l,
[0,64,128];
[0,128,255];
c l 7 ( l ,l ,
c l8 ( l,l,
[128.128.64];
c l9 ( l,l
[0,64,64];
c20(l,l
%*
% Assign numerical values for the first 10 coiore.
%
>^
fprinti(i>t,'tissues_properties.hin');
fprint{(l^t,Va!ues o f properties for tissuesW);
^rintf(ipt,"\n');
epsr = freespace(l);
sigma - freespace(2);
tirintfOJhWdefme FREE_SPACE_EPSR
%7.5ftn',apsr);
j|)rintff^t,'#define FREE_SPACE_SIGMA %7.5fm',sigma);
epsr = lungs(l);
sigma = lungs(2);
t>rintfl;!fpt;#define LUNGS_EPSR
%7.5fai',epsr);
tjr in tf ^t,'#def!ne LUNGS_SIGMA
%7.5f\n',sigraa);
epsr = b lood(i);
sigma = blood(2);
126
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
APPENDIX E
FDTD SIMULATION CODE
127
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
ihy == Res__Mem_T_hpcfloat(IE,JE);
Filename: 2dfdtdfepc03.c
2D TM simulation (Ez, Hx, Hy).
if(!only_Ezlmexwr)
{
amp___ez = Res_Mem_T_hpcfioat^E,JE);
#mctade <stdlib.h>
#mclude <stdio.!i>
#inciude <math.h>
#inc!yde <string.h>
arap = Res_Mera___T__bpcfloat(IE,JE);
p h i e = Res_Mem_T_hpcfioat(IE,JE);
real_pt ~ Res_Mem_T_hpcfloat(IE,JE);
imag___pt = Res_M8ra_T__hpcfloat(IE,JE);
ave__ampez = Res_Mem_T_fipcfloat(IE,JE);
#define NFREQS I
one_step_back_ez = Res_Mem_T_hpcfloat(IE,JE);
two__step_back__ez = Res_Mem_T_hpcfloat0E,JE);
typedef float T__hpcfloat; / ’^ Later use this to change T_hpcfloat to
double */
#include
#inc!ude
#include
#mclude
#include
#include
#include
#include
#include
#inc}ude
#include
#mclude
#incliide
peaks_coimter = Res_Mem_int(IE,JE);
}
"tissues.h"
"tissuesj>roperties.h”
"ernhpc.h"
'’dimensionshpc.h"
"mischpc.h"
"pmiiipc.h"
"fdtd2dhpc.h"
''save_2dfdtdl^c.h’'
"readjjararnhpc.h"
"initializehpc.h"
"ampEz__sshpc.h"
"mem_mgmt,h"
'ave_air^Ez.h"
pmpID = Res__Mem_int(IE,JE);
/* Initializing */
Initialize_Arrays3(ez_inc,hx_inc,ez,dz,hx,hy,ihx,ihy,iz);
Initialize__Arrays_ID(propID,constGa,constGb);
Rsset_AmpEz(amp_ez);
if(!only_Ezlinexwr)
{
Initialize_Arra>^_Fomier(realjDt,imagj[5t,amp,phase, //
real_in,imag_in,anp_in,phase_in,freq,arg);
}
ifi^!on1y__Ezlinexwr)
(
mainQ
/* Initializing Average Arraj^ */
Reset__fArray(ave_ampez);
Reset_fArray(one__step_back_ez);
Reset_fArray(two_step_back_ez);
Reset_iArray(peaks_counter);
{
T_hpcfloat **dz, **ez, *'*bx, '*'*hy, **iz, **ihx, ‘®'*ihy, **amp_ez;
T__hpcfloat gi2[IE], gl3[IE], fil[IE ], fl2[][E], fO[IE];
T_hpcfloat g}2[JE], gj3[JE], f]l[IE ], fj2[JE], §3[JE];
T_hpcfloat ez_m.c[JE], hx_mc[JEJ;
T_hpcfloat ez__inc_low_ml, ez_inc_low_m 2, ez__inc_high_mi,
sz_inc_high_m2;
int n, nsteps, T;
T_hpcfloat ddx, dt, tO, spread, source_lrequeiicy;
FILE *^ poin ti;
in tp o si_ x ,p o si_ 3 r;
int axis, initx, inity, endx, endy;
int movie, movie_period, steps, npml;
char dummy, filenamefD[80];
int T_ss, 1ast_period;
int **propID;
T_hpcfloatconstGa[NUM BER_OF_TISSUES],
constGbpWMBER_OF__TISSUES];
long hit size_axrays, size_dz;
}
/* Read parameters */
Read_Parameters(&ddx, &source_frequency, &tO, &spread,
&npml,
//
&movie, &movie_period, &nsteps, &axis, &initx, &inity,
&endx, &endy);
Read_ID__Filename(fiIenameID);
Read_Prop__j[D(propID,filenameID);
Read_Steady_State(&last_j>eriod);
Read__Average_Steady_State(&start_averaging_step);
Read_Only_Eziine(&only__Ez]inexwr);
Read_Start_Saving_Ezlmex(&stait_saving__Ezlmex);
/* Time step - Coiirant condition (Sullivan) */
dt = ddx/{2^ C );
T _ S ^ f i o a t ^ * r e a l j 3 t , * * im a g _ p t, * * a n ^ , * ’^ p h a s e ;
TJipcfloat rea!_in[NFREQS], imag_in[NFREQS],
3m p _ii^ F R E Q S ], phase_in[NFREQS];
T_hpcfioat foq[N FR EQ S], arg[NFREQSJ;
Calc_Constants(constGa,constGb,dt);
/* R eset 8 C for incident field */
ez_inc_low__m2 - 0.0;
ez_m c_!ow_m l = 0.0;
T__hpcfloat *^ave_arapez, **one_step_back_ez,
^*two_step_back_ez;
int *^eaks_coim ter, start_averaging__step;
int on!y_Ezlinexwr, start_saviDg_Eziraex;
FILE ^fpEziineswr,
ez_inc_higJL.^~
ez_inc_high_m l = 0.0;
Calc_PM L_Param eters(^2,gi3,f!l,fi2,fi3,g)2,gj3,fjl,§2,§3,npral);
Reserve menmry */
dz —Res_Mem_T_hpcfloat(IE,JE);
ez - Res_Mein_T_hpcfloat(iB,JE);
hx = Res_Mem_T_lijcfloat(IE,JE);
hy = Res__Mem_T__l^cfloat0E,JE);
iz = Res__Mem_T_h^float(IE,JE);
ihx = Res_Mem_T_hpcfloat(IE,JEX
/* Parameters for the Fourier Transforms
freq[0] = soiirce_frequency;
freq[l] = 300.e6;
fieq[2] = 700.e6;
for(m=0;n<NFREQS;n-H-)
128
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
if( T >= statt_saving__Ezlinex)
Save_Ez_linex_wr(fi5Ezlinexwr,ez,axi3,mitx,iiiity,endx,end.y);
i^(fppointi=fopen("Ezpolntr',"w"))-=NULL)
( printfi^'*\nThe file EzpointI could not be opened");
exit(0);
ifi^m ovie=YES)
{
}
i^steps==iBovie_period)
(
Save_Ez(ez);
steps = 0;
if(!onIy_E 2 linexwr)
{ Fourier_Amp_Phase_Inc(amp_iiii,phase_in,real_in 5imag_!n};
if((%Bzlinexwr=fopsn("Ezliiiexwr","w"))=NULL)
{ prmtfi^"\nl1ie file Eziinexwr could not be opened");
exit(0);
Save_Fourier_Ainp_Phase_Space(amp 5phase,ainp_in,phase_in,real_
pt,iinagj3t);
/* Ez points coordinates to save. */
p o sl_ x = IC;
p o sl_ y - inity;
Save Fourier Line(axis,amp.phase.amp in.pbase imreal pt.imag pt
.freq, //
}
}
/* Initialize time (iterations) ®/
T = 0;
steps = 0;
T__ss = nsteps - last_period;
mitx,inity,endx,endy);
for(m= 1;n<=nsteps;n+~^)
ifi^!oniy_Ezlinexwr)
{
{
T = T+1;
/* prmtf("%d\n‘',T); */
ifi; T > T_ss )
Update_AmpE 2 (ez,amp_ez);
/* — Start o f tbe Main F D ID loop — */
if( T >= start_averaging__step)
Caic_Inc_Ez(ez_mc,hx__inc);
/* printf("%d ”,T); */
Update_Ave_AmpEz{ez, ave__an ^ z, one_step_back_ez,
two_step_back_ez, //
peaks_counter);
if(ionly_Ezlinexwr)
{
Fourier_Inc_Bz(reaI_in,imag_in,arg,ez_inc,T);
}
if( T > = (start_averaging_step - 2 ) )
/* ABC for the incident buffer
ez_inc[0] = ez_incJow _m 2;
ez_inc_low_m 2 = ez_inc_low _m l;
ez_m c_iow _m l = ez_inc[lj;
ez_inc[JE-l] = ez_inc_high_ni2;
ez_inc_high_ni2 == ez_inc_high_ml;
ez_inc_high_m l = ez_inc[JE>2];
(
Copy_fArrays(one_step_back_ez 5 two_step_back_ez);
Copy_fArrays(ez, one_step_back_ez);
- End o f the main FDTD loop —
Calc__Dz(dz,gi3,gj3,gi2,gj2,hy,hx);
fclose(^5pointl);
fclose(^Ezlinexw r);
Save_Ez(ez);
Update_Source_Sinissoidal_PW(ez__inc,dz,dt,T,source_ffequency);
Update__Source__FW(e2 _mc,dz,dt,T, to,spread); */
/* Update_Line___Source_Pulse(dz,dt,T,tO,spread,lsx,lsy); */
/* Update_Line_Source_Sine(dz,dt,T,lsxdsy,source_&equency); *I
if(!on]y_Ezlinexwr)
{
Save_AmpEz(amp_ez);
Save__AmpEz_Line(axis,amp_ez,initx, inity,endx,endy);
Save_Ave_AmpEz(ave__ampez);
Save_Ave_AmpBz_Lme(axis,ave_ampez,initx,inity,endx,endy};
Fourier_Amp_Phase_Inc(amp_in,phase_m,realJn,miag_in);
Calc_Inc_Dz(dz,hxJnc);
Calc_Ez__ID(ez,dz,iz,propIDjConstGa,constGb);
iI(!only_Ezlinexwr)
Save_Fourier_Arap_Phase_Space(amp,phase,ainp_in,pliase_in,reai_
pt,imag_pt);
Fourier Ez(ez,real pt,imag pt.arg.T);
}
Save_Fourier__Line(axis,amp,phase,amp_m,phase_in,reai_pl,imagj3t
,freq, //
Caic_Inc_Hx(hx_inc,ez_inc);
C alc_H x(ez,hx,ihx,ezJnc,fil,fi2,fi3);
Calc_Hy(ez,hy,ihy,ez_inc,fi 1,l!2,fi3);
Save_Ez_Pomt(§3point 1,ez,pos l_x,p os l_y);
imtx,inity,ei5dx,endy);
Save__Fourier Amp_Phase_Inc(amp_m,phase_in,real_in,imag_in,fre
q );
*/
if(only_Eziinexwr)
129
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
#define
#define
#defm s
#defme
Sdefme
#defins
#define
#define
3ave_Fomer_Anp_Phase_Ax!s(am p,pbase,am p_in,phase_m ,realj5
t,imag_pt,freq,ga); V
Saye_Fourier_Xaxis(amp,phase,anip_in,phase_in,real_pt,iniag_pt,fi-e
q,ga);
}
BLOOD_EPSR
BLOOD_SIGMA
BONES_EPSR
BONES_SIGM A
HEART_WALL_EPSR
HEART_WALL_SIGMA
iNT_FLUID_EPSR
INT_FLUID_SIGMA
#define ESOPHAOUS_EPSR
printl('T = %d \n",T);
prmtl("Sirau!ated tsine = % i0 .9 f
naBosecs.\n",(T_tocfloat)(T*‘dt)* i e9);
printf("Siimiated FDTD space = % 10.9fiiix% 10.9fin
.Vi'',(T_hpcfioat)(IE*ddx),(T_hpcfloat)(JE'*'ddx));
printf("Simulated FDTD space = %d cells x %d cells .\n'',IE,JE);
printfC'Size o f T_hpcfioat is %d\n",sizeofi^T_hpcfloat));
64.79700
#defme ESOPHAGUS_SIGM A
1.23160
#defilie INT_MUSCLES_EPSR
54.81100
#define INT_MUSCLES_S1GMA
0.97819
#defm e THORACIC_WALL_EPSR 28.74150
#defme THORACIC__WALL_SIGMA 0.50800
#define TUMOR_EPSR
56.51340
Sdeiine TUMOR_SIGMA
1.29111
printf('TE = %d, JE = %d, IE*JE = %d \n*',IE^E,IE"‘JE);
size__dz == sizeofiTJbpcfloat)’^IE*JE;
printf("The size o f dz is %d\n",size_dz);
piintf("Since w e have 9 dz-iike arrays, then memory req. =
%d^n'’,9*size_dz);
em.h
EM constants.
/* Free memory
/*Free space permittivity *!
Free_Mem__T_hpcfloat(d2 jE,JE);
#define EPSZ 8.854l9e-12
Free_Mem_T_hpcfloat(ez,IE,JE);
#defm e C 3e8
Ftee_Mem_T_hpcfloat(hx,IE,JE);
Free_Mem_T_hpcfloat(hy,IE,JE);
Free_Memjr_hpcflo8t(iz,IE,JE);
Free_Mem__T_hpcfloat(ihx,IE,JE);
Free_Mera_T_l^cfloat(ihy,IE,IE);
Free__Mem_in^ropID,IE,JE);
dimensionshpc.h
These dimensions correspond to a
1-GHz analysis
/♦ Dimensions o f FDTD space */
#defineIE 158
#defme JE 158
/* TotaPscattered field boundaries */
# d efin e lA 23
#define IB (IE-IA-1)
if(!on!y_Eziinexwr)
I
Free_Mem_T__hpcfioat(an^_ez,IEJE);
Free_Mem_T_hpcfloat(ave__an:^ez,IE,JE);
Free_Mem_T_hpcfloat(one_step_back_ez,IE,JE);
#defineJA23
#defineJB(JE-JA-l)
Free__Mem_T_hpcfloat(two_step_back_ez,IE,IE);
Free_Mem_int(j>eaks_counter,IE,JE);
}
}
/* Coordinates ofF D T D space center */
#defme IC IE/2
#define JC JE/2
tissues.h
Identifiers for tissues
#define NUMBER_OF_TISSUES
#defme
#defme
#define
#define
#de& ie
#defme
#define
5 i .06500
1.58290
16.60784
0.26935
59.29000
1.28360
68.87500
1.66730
FREE SPACE
LUNGS
BLOOD
BONES
HEART_WALL
INT_FLUID
misc.h
10
#de{m ePI
//define TRUE
Sdefine FALSE
#defineY E S
S defineN O
0
1
2
3
4
5
ESOPHAGUS
#define INT_MUSCLES
6
7
#define THORACIC_WALL
#define TUMOR
8
9
3.14159
1
0
i
0
void Beep(void);
void
B e^ (v o id )
{
prii!t{("\a");
}
tissues_properties.h
Values o f properties for tissues
These vatoes are for 143Hz analysis
# d e fc e FREE_SPACE_EPSR
t.OOOOO
#defme KREE_SPACE_SIGMA
Sdefme LUNGS_EPSR
#define LUNGS SIGMA
0.00000
21.82500
void Calc_PML_P 2 rameters(T__hpcfioat gi2[TE], T Jipcfloat gi3[IE],
//
T_hpcfloat filpE], T_hpcfloat fi2[IE], T_hpcfloat f!3[IE], //
0.47406
130
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
T^Jipcfloat gj2[)E]j T_hpcfloat gj3[JE], T__bpcfloat
T _l^ ciloat §2[JE], T_h|3cfloat §3[JE], int npml);
ljl[3E-2-j] = xn;
g if j ] = !.0/(l.(M-xr.);
P [J E - 2 - j ] = 1.0/(i.0+xn);
P [ j ] ‘= (i-0 -s n )''(!-0 + xn);
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn):
//
1
void
Ca!o_PML_Paiamete!^T_hpK:float gi2[lE],TJipcfloat gi3[IE],
T_hpcfloat fil [IE], T_hpcfloat fi2[IE],T_hpcfloat Ss3[II3], //
TJipcfloat gj2[JE], TJhpofloat gj3[JE],T_hpcfloat §t[JE ], //
T_hpcfloat 52[JE], T_hpcfloat fj3[JE], int npml)
$ if:^li;^.It$$ $^
$tS*-4::$s -i:^ ^
^ ^ =4: tf:4s 1#$ ^
if:$ ^
^
- Filename: fdtd2d.h
{
int ij;
/* int npml; */
T_hpcfloat xn, xxn,
V
/* FDTD Analysis */
void Calc_Inc___Ez(T_hpcaoat ez_inc[JEj,T_hpefioat hx_mc[IE]);
void Calc_Dz(T_hpcfloat
T_hpcfloat gi3[IE],T_hpcfloat
xnutn, xd;
m m ji
for(i=0;i<IE;i++)
{
T_hpcfloat gi2[IE],T__bpcfloat gj2[JE],T_hpcfloat **hy, //
T_hpc£loat *’*'hx);
gi2[i] = !.0;
gO[i] = 1.0;
void Calc_Inc_Dz(T_i^cfloat **dz,T_hpcfloat hx_inc[JE]);
void Calc__Ez(T_hpcfloat ez[IE][JE],T__hpcfloat dz[IE][JE}, //
T_hpcfloat iz[IE][JE],T_hpcfloat ga[IE][JE], //
fil[i] = 0.0;
fi2[i]= 1.0;
fi3[i]= 1.0;
T_hpcfloat gb[IE][JE|);
void Calc_Ez_ID(T_hpcfloat **ez, T__hpc£ioat
//
T__hpcfloat **iz, int **^ropID,
11
T__i^cfloat constGapSrUMBER_OF_HSSUES], //
T__hpcfloat constGb[NUMBER_OF_TISSUES]);
void Calc_Inc_Hx{T__hpcfloat hx_inc[IE],T_hpcfioat ez__mc[JE]);
void Caic__Hx(T_hpcf!oat
T_hpcfloat *‘*‘hx,
//
T_hpcfloat **shx,T_bpcfloat ez_mc[JE], II
T_hpcfloat fiip E ], T_hi>ct!oat ^2[JE] J „h p cfloat g3[JE]);
void Calc_Hy(T_hpcfloat **ez, T_hpcfloat **hy,
//
T _l:^ float **ihy, T_hpcfloat ez_inc[JE], //
T_hpcfloat f 1 [JE], T_hpcfloat fi2[IE], T_hpcfloat fi3[IE]);
)
for(j“ Oa<IE;j++)
{
gj2[j] = I.O:
g)30] = 1.0
giO ] = o.o
P M = 1.0
P D ] = i-O;
}
/* printf("Number o f PML cells —> ”); *!
I* scanf("%d", &npml);
*/
/* Source ^5datmg */
void Update__Source_PW(T_hpcfloat ez_mc[JE],T_hpcfloat
dz[IE][JE], //
T_hpcfIoat dt, int T,T_hpcfloat tO,T_hpcfloat spread);
voidUpdate_Source_Sinusoidal_PW (T_hpcfloat
ez_inc[JE],T_hpcfloat **dz, //
T___hpcfloat dtjint T,T_hpcllo3t freqi^ncy);
voidUpdate_Lir^_Source(T_hpcfloat dz[IE][JE],T_hpcSoat dtjint
for(i=0;i<=npml;i-H-)
{
xnum = npml - i;
xd = npml;
xxn = xmim/xd;
sn = 0.33 * pow(xxn,3.0);
gi2[i] = 1.0/(1.0+xn);
gi2[IE-l-i] - 1.0/(1.0Txn);
gl3[i] = (1.0 - xn)/(1.0 + xn);
gi3[IE-i-l] = (1.0-xn)/(1.0+xn);
xxn = (xmim - 0.5)/xd;
xn = 0.25 * pow(xxn,3.0);
f il [ i| = xn;
fil [IE-2-i] = xn;
fi2[i] = 1.0/(1.0+xn);
fi2[IE -2-ii= 1.0/(l.(H-xn);
fi3[i] = (I.O - xn)/(1.0 + xn);
fl3[IE-2-i] = (1.0 - xn) / (1.0 + xn);
IJl
TJhpofloat tOjTjhpcfloat spread,int IsXsint Isy);
void UpdatejLiiie_SourcejPulse(TJipcfloat dz[IE][JE],T_hpcfloat
dt.intT,//
Tjhpcfloat tO,T__hpcfloat spread,int !sx,int Isy);
void Update_Line_Soiirce_Sine(T_hpcfloat dz[IE][JE],T__hpcfloat
dt,int Tjint Isx, //
int lsy,T_hpcfloat freqi^ncy);
Fourier Analysis'*'/
void Fo'uiier_Inc_Ez(T_hpcf!oat real_m[NFREQS],Tjhpcfloat
imagjin[NFREQS], //
Tjhpcfloat arg|NFREQS],T_hpcfloat ez_inc[JE],int T);
void FourierjEz(T_hpcfioat **ez,T__hpcfloat
//
Tjhpcfloat **iim gjst,T _hpcf!oat arg[NFREQS],int T);
void Fourier_Airyp_Pliase__Inc(TJ^cfloat
an^_in[NFR EQ S],T_l^cfloat p h a se_ ii^ F R E Q S ], //
Tjhpcfloat real_in[KFREQS],T_iq5cfloat
imag_m[NFREQS]);
void Fourier_Lme__Source(Tj!ipcfloat real_m[NFREQS],TjF^cfloat
imag^m[NFEEQS], //
T_i^cfloat arg[NFREQS], T_hpcfloat ez[iE][JE],mt T, //
int Isx, int Isy);
void Ca!c__Constants(T_hpcfloat
}
for(j=Oy<=npmly-H-)
{
xnum = npml -j;
xd - npml;
xxn = xmim/xd;
xn = 0.33 * pow(xxn,3.0);
g j2 0 ]= 1.0/(!.0+xn);
gj2[JE-l-j] = !.0/(1.0+xn);
gj3Q] = (1.0 - xn)/(1.0 + xn);
gj3[JE-j-I] = (i.O - xn) / (1.0 +
xxn = (xnum - 0.5)/xd;
xn = 0.25 * pow(xxn,3.0);
xn);
constGa[NUMBER_OFjTISSUES], //
Tjhpcfloat constGb[NUMBER_OF_TISSUES], T_hpcfloat dt);
giO] = xn;
131
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Calo_Inc_Dz(T_hpcfloat ®'*dz,T_!ipcf!Qat hx_it!c[JE])
V
i n t i;
Calc_!nc__Ez(T_iipcfloEt ez__inc[JE], T__hpcfioat hx__inc[JE])
for(i=IA;i<=IB;i++)
{
{
dz[i][JA] = dz[i][JA] + 0.5 * hxJtic[JA -l];
dz[i][JB] = dz[i][jB] - 0.5 * hx_incpB];
intj;
for(i=iy<JE;j-H-)
{
€z__inc{j] = ez_inc|j] + 0.5 * (hx_inc[j-l] - hx_inc[j]);
}
}
Calc_Ez(T_hpcfloat ez[IE][JE],T__lipcfioat dz[IE][JE],T__hpcfloat
iz[IE][JE], T_hpcfloat ga[IE][JE],T hpcfloat gb[IE][JE])
{
void
Fomier_Inc_Ez(T__hpcfloat reaHnfNFREQ S], T_l^cfloat
imag_iii[NFREQS], //
T hpcfloat arg[NFREQS], T_hpcfioat ez inc[JE], int T)
int iJ;
/®Infd2d_3.3.c
{
foi(j=Oa<!Ey++)
{
intm;
for(i=0;KIE;i++)
for(m.=0;m<NFREQS;m-H-)
{
for(i=ly<JE -ly-H -)
real_in[m] = real_in[m] + cos(arg[m]*T) * ez_m c[JA-l];
imag_in[m] = imag_in[m] - sin{arg[m]*1} * ez_iuc[JA-l];
{
for(i= 1;i<IE-1; i++)
{
ez[i]0]=ga[i]D]Vdz[i][)]-iz[i]U]);
}
}
M i]D] = iz[i]D] + g M i]0 ]* e z [i]0 ];
}
void
Calc_Dz(T_hpcfloat *‘*'dz, T__hpcfloat giSpE], T___hpcfloat gj3[JE], II
T_hpcfioat gi2[IE], T_hpcfloat gj2[JE], T_hpcfloat *’*'hy, //
T_i^cfloat **Ik)
/* Set the Ez edges to 0, as part o f the PML “I
{
for(j=Oy<JE-ly++)
{
int iJ;
ez[0][j] = 0.0;
ez[IE -i][j] = 0.0;
for(j=ly<JEy+4-)
/* In Sullivan is IE! */
(
for(i-l;i<IE;i-i-i-)
{
foifi=0;i<IE-I;i++)
{
dz[i][j] = gi3[i] * gj30] * dz[i]D] + gi2[i] »gj20] * //
0.5 * (hy[i]0] - hy[i-l]D ] - ta[i][j] + hx[i]0-l]);
ez[i][0] = 0.0;
ez[i][JE -i] = 0.0;
}
}
}
void
Update___Source_FW(T__hpcfloat ez_!nc[JE],T_hpcfloat
dz[IE][JE],T_hpcfloat dt,mt T, //
T_hpcfloat tO,T_hpcf1o2t spread)
Fomier_Ez(T__hpcIloat ^*ez,T__hpcfloat **real_pt,
//
TJipcfloat **imag_pt,T__hpcfloat arg[NFREQS],int T)
{
int ij,m ;
{
T_l^cfloat source;
f o r ( i= 0 ;|< J E ;j- H - )
{
source = sm (2^PP700=*^le6W T);’^/
source = exp(-0.5*(pow((tO -T )/sprei,2.0)));
for(i=0;i<IE;i++)
/* In Sullivan is JE
{
for(m=0;in<NFREQS;m4H-)
{
source; */
ez_inc[JA-4] = source;
/* e2_inc[3] =
real_pt[i]0] “ real_pt[i][j] + cos(arg[m] ®T) * ez[i]fl];
inrag_pt[i]|j] = imag_pt[i][j] - sin(aig[m] » T) ® ez[i]Q];
/* According to fdld_2.2.c this should be hut in fd2d 3.4.c is + */
void
132
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
f o r(i= J A -j< = J B ;j+ + )
{
}
by[!A -!]0] = h y P A -l]0 ] - 0.5 » e z jn c[j];
by[IB][ij = hy[IB]0] + 0.5 » ezJncO ];
}
}
void
Calc_Inc_Hx(T_hpcfloat hx_inc[JE], T_bpcf!oat ez_inc[JE])
Calculate the Fourier amplitiide and phase o f the incident pulse
j;
void
Fourier__Afs:^_Pliase__lGc(T_hpcfioatanip_in[N FREQS],T_hpcfloat
phase_m[NFREQS], //
T_hpcfloat real_m[NFREQS]jT_hpcfloat
imag_in[NFREQS])
for(i=Ou<sEa-*^)
!
hx_inc[j] = hx_inc[j] + 0.5 * (ez_inc[i] - ez_m c|j+ l]);
for(m=0;m<NFREQS;m-f-+)
{
amp_m[m] = sqrt(pow(real_in[m],2.0) + pow(imag_in[m],2.0));
phase_in[mj = atan2(imag_in[m],real_in[m]);
/* prmt£(''%d Input Pulse : % S.4f % 8.4f % 8.4f
%7.2fin",m,real_in[m], //
void
Caic_Kx(T_hpcfloat *^ez, T_hpcfloat **hx, T_hpcfloat ’^^ihx, !!
T_hpcfloat ez_mc[JE], T_hpcfioat fil[IE ], T_hpcfloat 02[JEJ,
//
T_hpcfloat 53 [JE])
}
}
{
in tij;
T_hpcfloat c«rl_e;
imag_in[m],ainp_m[m],( 180.0/PI) * phase_in[m]); •/
for(j=05<JE -1:j++)
{
for(i=0;i<IE;i-H-)
{
Fourier_Line_Source(T_hpcfloat real_in[NFREQS],T_hpcfloat
imag_inp4FREQS], //
T_hpcfloat arg[NFREQS],T_lipcfloat ez[IE][JE],int T,int
lsx,int Isy)
curl_e = ez[i][j] - ez[i]0+l];
ihx[i][j] = ihx[i][j] + flip] * curl_e;
hx[i]0] = §301 * hxplOl + § 2 0 ] * 0.5 *(curl_e + ihx[i]0]);
for(m=0;m<NFREQS;m++)
{
/* Incident Hx values *1
rea!_in[m] =
real_in[m] + cos(arg[m]*T) * ez[lsx][lsyj;
iraag^in[m] = imag_m[m] - sin(arg[m ]*r) * ez[lsx][lsy];
for(i=IA;i<=IB;i++)
{
hx[i][JA-l] = hx[l][JA-l] + 0.5 » ez_mc[IA];
hx[i][JB] - hx[i][JB] - 0.5 ®ez__inc[JB];
void
Update__Line_Source(T_i3pcfloat dz[iE][JE],T_hpcfloat dt,int
T,T_hpcfioat to, //
T_hpcfloat spread,int lsx,int isy)
1
{
Calo_Hy(T_hpcfloat **ez, T_hpcfloat *®hy, T_hpcfloat *®ihy, //
T_hpcfloat ez_mc[JE], T_bpcfloa! 5 1[JE], T_hpcfloat fi2[IE], II
T_hpcfloat f6[IE ])
T_hpcfloat source;
source = sin(2*PI*700* le6Mt*T);
{
/* source = exp(-0.5*(pow((t0-T)/spread,2.0))); V
dz[isx][isy] = source;
in tij;
T_hpcfloat curl_e;
}
for(j=Oy<=JE-ly++)
{
for(i=0;i<IE- i ;i++)
{
void
Update_Source_Sinusoidal_PW(T_lipcfloat ez_inc[JE], TJipcfloat
curl_e = ez[i+l][j] - ezpjO];
ihy[i][j] = ihy[i][j] + §10] * ourl_e;
hyPlO] - fi3[i] * hyplOl + fi2[i] * 0.5 * (curl_e + ihy[i][j]);
T__hpcfloat dt, int T, T_hpcfloat frequency)
}
}
T_bpcfloat
source, w ,
t, wt;
t-T * d t;
w = 2 * PI * frequency;
/ • Incident Hy values */
133
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
w t - w * t;
source = sin(wt);
ez_inc[JA-4] = source;
for(i=0;i<IE-S ;i++)
{
}
ez[i][0] = 0.0;
ez[i][JE -l] = 0.0;
Update_Lme_Soun:e_Pulse(T_hpcfloat dz[IE][JE],T_hpcfloat dt,int
T,T__hpcfloat to, //
T__]^float spread,int isx,int Isy)
void
Ca!c_Constants(T_hpcfloat coiistGa[NUMBER__OF_TISSUES], //
T_hpcfloat const£jb[NUMBER_OF_TISSUES],
T_l:^}cfloat dt)
{
T_hpcfloat source;
source - exp(-0.5*(pow{(t0-T)/spread,2.0)));
dz[lsx][lsy] = source;
!
constGa[FREE_SPACE] ~ i.O/(FREE SFACE_EPSR +
(FREE_SPACE_SIGMA*dt/EPSZ) );
constGb[FREE_SPACE] = FREE_SPACE_SIGMA»dt®PSZ;
constGa[LUNGS] = 1.0/(L im G S^E F SR +
(LUNGS_SIGMA*dtffiPSZ));
constGb[LUNGS] = LUNGS_SIGMA*dt/EPSZ;
Update_Line_Soiirce_Sme(T_hpcfioat dz[IE][JE],T_l^cfloat dt,int
T,int IsXjint Isy, T_l:q?cfloat frequency)
oo 0 StGa[BLOOD] = 1.0/(BLOOD_EPSR +
(BLOOD_SIGMA*dt/EPSZ) );
constGb[BLOOD] = BLOOD_SIGMA»dt/EPSZ;
constOa[BONES] = 1.0/(BONES_BPSR +
{
T__hpcfloat source;
source = sin(2*Pl*frequency*dt*T);
dz[lsx][lsy] = source;
(BONES_SIGMA*dt/EPSZ) );
}
constGb[BONES] = BONES_SIGMA*dt/EPSZ;
constGa[HEART_WALL] = 1.0/(HEART_W ALL_EPSR +
(HEART_W ALL_SIGMA*dt/EPSZ));
con stG b piE A R T _W A Ii] = HEART_WALL_SIGMA*dtffiPSZ;
constOa[INT_FLUID] = 1.0/(INT_FLUID_EPSR +
(INT_FLUlD_SIGMA*dt/EPSZ) );
constOb[INT_FLUID] = INT_FLUID_SIGMA*dt/EPSZ;
constGa[ESOPHAGUS] = 1.0/(ESOPHAGUS_EPSR +
(ESOPHAGUS_SIGM A«dt/EPSZ));
void
C3lc_Ez_ID(T_hpcfloat *®ez, T_hpcfloat **dz, T_hpcfloat ®*iz, //
int **propID, T_hpcfloat
constGa[NUMBER_OF_TISSUES], //
T__hpcfloat confitGb[NUMBEK_OF_TISSUES])
i
constGb[ESOPHAGUS] = ESOPHAGUS_SIGMA*dt/EPSZ;
constGa[INT_MUSCLBS] = 1.0/(INT_MUSCLES_EPSR +
(INT_MUSCLES_SIGMA*dt/EPSZ) );
constGb[INT_MUSCLBS] = INT_MUSCLES SIGMA*dt/EPSZ;
int ij.id;
T Jipcfloat value_ja, value _^b;
P In f<J2d_3.3.c
fot(j=0;j<JE;j++)
C0MtGa[TH0RACIC_WALL] = 1.0/{THORACIC_WALL_EPSR
+ (THORACIC_WALL_SIGMA*dt/EPSZ));
constGb[THORACIC_WALLJ =
«/
THORACIC_WALL_SIGMA*dt/EPSZ;
con3tGa[TUMOR] = 1.0/(TUMOR_EPSR +
for(i=0;i<IE;i++)
(TUMOR_SIGMA»dt/EPSZ) );
constGb[TUMOR] - TUMOR_SIGMA»dtffiPSZ;
/* g a = S.O/(epsr + (sigina*dt/EPSZ) );
gb = sigma*dt/EPSZ; */
foi(j= l;j<JE -la++)
{
for(i= 1;i<IE-1; i-H-)
}
id = propID[i]0];
printlj"i = % d, j = %d, id = %d‘in",ij,id );
value_ga = coastGapd];
printi{''value_ga = %Pn",value_ga);
value_gb = constGbpd];
/*
prmtf("value_^ = %f\n",va!ue__gb);
/*
*/
save_2dfdtd.il
*/
/* #defme
Save , Fourier Yaxis(amp,phase,amp in,phase in,real pt,imag ptfre
q ,g a )//
*/
ez[i][j] = vai«e_ga * (dz[i][j] - iz[i]0]);
iz[i]0] = iz[i]D] + valu£_gb * ez[i]0]:
Save_Fourier_An^_Phase_Axis(axnp,phase,anip_in,phase_in,real_p
t ,//
iinag_pt,freq,ga) */
}
/* Set the Ez edges to 0, as part o f the PML */
void Save_Ez(T_hpcfloat **ez);
void Saye_Fourier_An:^_Phase__iiic(T_hpclloat ani|5_m[NFREQS],
T_hpcfloat p h a se_ in [N F ^ Q S ],//
T_hpcfk>at real_in[NFREQS], T_hpcfIoat Lmag_m[NfREQSJ,
T_l:^cfloat fl-eq[NFREQS]);
void Save_Fourier_Amp_Phase_Space(T_i:^cEoat **amp,
T_hpcfloat **phase, //
fo r(i= O y < J E -iy + + )
{
ez[0](j] = 0.0;
ez[IE -l]0] = O.O;
134
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
T_bpcfloat anipJn[NFREQS], TJipcfloat
phasejBp^FR EQS],
fl
T_hpcfloat **i’eal_j)t, T_bpcfloat *^*imagjpt);
void Save_Fouri€r_Amp_Phase_Axis(T_hpcfloat amp[IE][JE],
T_hpcfloat pbase[IE][JE],
//
T_hpcfloat amp_in[NFREQS], T_hpcfloat
phaseJn[N FR EQ S],
//
T__hpcfloat real_j5t[NFREQS][fE][JE], T_bpcfloat
im a^ t[N F K E Q S][IE ][JE ],
//
T ^ l^ flo a t freq^FR E Q Sj, T_hpcfloat ga[IE][JE]);
void Save_Ez_Point(FILE *^point, T_hpcfloat **ez, int pos_x, int
p o sj);
void Save_Ez_Y_axis(T_hpcfloat ez[IE][JE], intpos_x);
void Save_Foisrier_Xaxis(T_hpcfloat amp[IE][JE], T_hpc£loat
phase[IE][JE],
//
T^hpcfloat ampJn[NFREQS], T_hpcfloat
phase_m[NFREQS],
//
T Jipcfloat real_j)t[NFREQS][IE][JE], T_hpcfIoat
imag_pt[NFREQS][IE][JE]»
//
T_hpcfloat feq[N FR E Q S], T_hpcfloat ga[IE][JE]);
void Save_Fourier_Lme(int axis, T_hpcfloat
T_hpcfloat
'**phase,
void
S ave_Fourier_A m p_PhaseJnc(T J^float
amp___in[NFREQS],T_hpcfloat phaseJn[NFREQS], //
T__fapcfloat reai_m[N'FREQS],T__h|)cfloat
imag^m[NFREQS|,T_hpcfloat fi-eqJJJFREQS])
{
FILE
intm;
^ = fopen("fapip'V’w");
for(nv=0;m<NFREQS;m-H-)
{
amp_in[m] = sqrt(pow (realjn[m ],2.0) + pow(imag_in[m],2.0));
phase__jn[m] = atan2(imag_in[m],real_m[m]);
j^rintf^lpf’Frequency %d = % 8.2f MHz Real part = % 8.7f
Imaginaiy
part = % 8.7f Anplitude = % 8.7f Phase - % 8.7f \n",m,
&eq[m3*i.e-6 ,reai_in[ra],imag_in[m], am3J_in[m],
(180.0/PI) * phase_inCm]);
//
T_hpcfloat air^_in[NFREQS], T_hpcfloat
phase_in[NFREQS],
//
T_hpcfloat **real_pt, T_hpcfloat **in»g_pt,
T_bpcfloat freq[NFREQS], int initx, int inity, int endx, int
endy);
void Save_Ez__linex_wr(FILE *^Ezlinexwr, T_hpcfloat **ez, int
axis, //
int initx, int inity, int endx, int endy);
*/
}
fclose(^ );
}
II
Save the Fourier amplitude and phase o f the total field in the whole
FDTD
space. The phase is saved in file in degrees,
void
Save_Fourier__Amp_Phase_Space(T_hpcfioat **anp,T_hpcfloat
Write the E field out to a file *'Ez"
’^^phase, //
void
Save Ez(T_hpcfloat **ez)
T_hpcfloat amp__in[NFREQS],T_hpcfloat phase_in[NFREQS], //
T_hpcfloat **real_pt,T_hpcfloat **imag pt)
I
char fi!ename[20];
FILE
FILE
’^fpphase;
int ij,m ;
T_hpcfloat phase_deg, phase_rad;
r
printf("Enter filename to store Ez data:");
fflusli(stdin);
gets(filename);
foF(m=0;nKNFREQS;m-i-+)
{i^jiy==0)
{ ^ = fopenC'tan^i'V’w");
^ h a s e = fopenC'tphasel'Jw");
if((^ fo p en (fa eiia m e" w " ))= N U L L }
{ printi^"\nThe fiie %s could not be opened",filenarae);
exit(0);
}
sise !f(n F = i)
{ ^ - fopen("tarap2'Vw");
^ h a s e = fopen("tphase2","w");
}
if((^fopen("E zalpha","w "))=N U L L )
{ printf("\nThe file Ezaipha could not be opened");
exit(0);
eise iHm==2)
{ = fopen("tanip3","w'');
fpphase = fopen("tphase3'V'w");
}
for(j=Oy<JEa++)
forG==Oy<JBy-H-)
{
for(i=0;i<IE;i++)
{
fprinti{fp,'' % !5.14f',ez[i][j]);
}
W);
}
fclose{^);
)
{
for(i=0;i<IE;i++)
{
amp[i]jj] = (1.0/an5p_in[m]) * sqrt{pow{real_pt[i]0],2.C) + //
pow (inm gj)t[i][|],2.0));
55rintft),”% i0.9f ”,aiHp[i]0]);
phase[i][j] - ataa2(imag_pt[i][j].reai_pt[i]0]);
phasejrad = phase[i][j];
Save the amplitude and phase o f the incident pulse at frequencies o f
interest.
/* Show the phase values in range 0 to -360 "*/
135
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
I*
Converting phase t
phsse_deg = phase_rad*(i 80/PI);
printf("%2d %9.8f\B",JC-J,phase_deg);
fprintfl;fpphase," % 9.Sf \n",phase_deg); */
fprintf{(pphase," % 10 .9 f \n",phase_deg);
'<f(phase[i]0]>0-0)
{
phase_rad = pSiase[iJ[j] - 2*PI;
}
*r
/* Converting phase to degrees. */
phase_deg - phase_rad * (ISO,'PI);
/* l^mtf(% phase ”% 9.8f ",phase_deg);
%rintii;^phaae,"% 10 .9 f ",phase_deg);
fciose(fp);
fclosc(f^phase);
}
}
}
fprintfi[ip,” W );
fj)rintf(^phase," \n");
}
fclose(ip);
fclose(^phase);
void
Save_Ez_Point(FILE ’^'^point, T_hpcfloat
pos_y)
}
}
int pos_x, int
{
T_hpcfloat ez_value;
ez_yalue = ez[pos__x][pos_3 '];
^rmtf(^poiat,''% 15.14f\n", ez_value);
Calculate and save the Fourier amplitude and phase o f the total field
aioiig the y axis inside the stnioture.
}
/»:
void
Save_Fourier_Ainp_Phase_Axis(T_hpefloatamp[IE]PE],T_hpcfloat
phase[IE][JE], //
T_hpcfloat atnp_tnP'IFREQS],T_hpcfloat phase_in[NFREQS], //
T_hpcfloat real_pt[NFREQS][IE][JE],T_hpcfloat
imagj)t[NFREQS][IE][JE], //
T_hpcfloat freq|WFREQS],T_hpcfloat gapE][JE])
void
Save__E2 __Y__axis(T_hpcf!oat ez[!E][JE], m tpos_x)
{
«
FILE *fp;
int j;
T_hpofloat ez_value;
int jjitt;
FILE
*fpphase;
T_hpcfloat phase_rad, phase_deg;
fp = fopen("Ez_y_axis'',"w”);
foi(j=Oa<tE;j++)
for(m=0;tn<NFREQS;m++)
{
{
if (m = 0 )
{ ^ = fopenC'amply'Vw");
ij^ hase = fopen("phasely","w");
ez_value = ez[posjs][j];
fprintfi^^" %10.9fvi",ez_value);
}
}
}
else i f ( n i = l )
{ = fopen(”amp2y'',''w'');
Ipphase = fopen("phase2y","w");
Calculate and save the Fourier amplitude and phase o f the total field
along the X axis inside the structure.
}
else if(m— 2)
{ fp = fopen("ainp3y'V'w");
fpphase = fopen("phase3y","w'');
Save_Fourier_Xajds(T_hpcfioat amp[IE][JE ]5 T_hpcfioat
phase[iE][JE],
//
T_hpcfloat aH^_m[NFREQS]j T__hpcfloat phase__mp4FREQS], //
T J ^ f lo a t reai_pt[NFREQS][IEf[JE], T _ l^ flo a t
imag_pt[NFREQS][IEj[JE], H
T__hpcfloat &eq[NF^EQS], T_hpcfioat ga[IE][JEl)
}
primC'% 2d% 7.2fM Hzln",m,frsq[m]*l.e-6);
for(j=JAy<=JBy++)
{
int i,m;
FILE
“^^phase;
T_hpcfloat phase_rad, pbase_deg;
ii(ga[IC]0] < 1.00) /* Saving only the data within the structure. */
{
/* Outside the structure Ga is 1 (free space) */
atnppCJO] = (1.0/ainp_in[iii]) * sqrt(pow(rsal_pt[m][IC][j],2.0)
+ //
pow(iinagj)t[m][IC][jJ,2.0));
printfr% 2d% 9.SfW yC -j,anip[IC ]0]);
/* |jrmtf[fp," % 9.8f W ’,amp[IC]0]);
^ rintffts,” % 10.9fW ,amp[ICjO]);
/ • phase[IC][j] = atan2(imaz ptfmlflCiriLreal p tfm in ciril) - //
phase_in[m];
*/
phase[IC]0] = atan2(itmg_pt[m][IC]0],rea!_pt[ra]pC]0]);
/ • Save the phase values in range 0 to -360 */
phase_rad = phase[IC][j];
i% hase[IC][j]>0.0)
for(iiF=0;m<NFREQS;m+'t-)
{
if(n F = 0)
{ = fopen("amp!x'V'w");
^ p h ase = fopenC'phaselK’V'w");
}
else ifi^ m = i)
{ ^ = fopen(''an^2x'V'w");
^ p h ase = fopen("phase2x'V'w”);
}
{
phase_rad = phaseJICJO] - 2*PI;
S
else if(m— 2)
fopen("ainp3x",”w");
fpphase = fopen("phase3x"V’w");
136
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
printf("%2d % 7.2fM H zW 5m ,freq[m ]*l.e-6);
for(i=IA;i<=IB;i'H')
{
if(ga[i][JC| < 1.00)
amp[i][inity] - (1.0/amp_in[m])
sqrt(pow(realjpt[i3[mity],2.0)
//
pow(!m agj)t[i][im ty],2.0)};
/* ^rintf(^," % 9.8f V ,am p[i][im ty]);
% 10.9f \n",an:^[i][mity]);
/* phase[i][init>'] = atan2(imag_j3t[i][inity],real_pt[i][imty]) - //
phase_in[m];
phase[i3[iiiity] = atan2(imag_ptFi][inity],realj3t[i][inity]);
Saving only the data within the structure.
r-^ Outside the structure Ga is i (free s
amp[i][JC] ~ (1.0/am p jn [m j) * sqrt(pow(real_pt[m][i][JC],2.0)
pow(!mag_pt[m][i][JC],2.0));
printl("%2d%9.8f\ii",IC-Unip[i][JC]);
/*
% 9.8f \n",amp[i][JC]);
^ n n tf(ip ," % i0.9fV ,am p [i][JC ]);
/* phase[i][JC] = atan2(imagj5t[m][i][JC],reaI__pt[m][i][JC]) - //
phase_in[m];
*/
phase[i][JC] ~ atao2(imag_pt[ra][iJ[JC],reai_pt[m][i][JC]);
{
/* Save the phase values in range 0 to -360
phase_rad ==phase [i] [inity];
lf(phase[!][mity]>0.0)
{
phase__rad = phase[i][jnity] - 2*P1;
/* Save the p h ^ e values in range 0 to -360 */
phase__rad = phase[i][JC];
if(phase[i][ICl>0.0)
Converting phase to degrees. */
p h ^ e_ d eg = phase__rad*( 180/PI);
/* lpniitf(i|)phase,'' % 9.8f \n",phase_deg);
fprintf(^phase," % 10.9f\n",phase_deg);
(
phase_rad = phase[i][JC] - 2*PI;
}
fclose(^ );
fclose(^phase);
} /* End o f FOR *!
/* Converting phase to degrees. '*'/
phase_deg = phase_rad*(l80/FI);
printfl;''%2d % 9.8f \n",iC-i,phase_deg);
/* fprintf(:^phase," % 9.8f \n",phase_deg);
^ r iiit^ ^ h a s e ,” % 10.9fV ,ph ase_d eg);
}
if(axis==2) /* A long Y */
{
for(m=0;m<NFREQS;m-H-)
{
}
fclose(^ );
fclose(^phase);
il( m = 0 )
{ <p = fopen(''amp_lmely'V'w");
4)phase = fopen(''phase_linely'’,"w'');
}
Calculate and save the Fourier amplitude and phase o f the Ez field
along a line, which can be either in X or Y axis,
axis: 1->X, 2->Y
else i f ( m = l )
( flj = fopen("amp_line2y",''w'');
fpphase = fopen{"pliase_lme2y'',''w");
}
void Save_Fourier_Lme(int axis, T_hpcfloat **an^,
//
T_hpcfloat **phase, T_hpcfloat anq5_in[NFREQS],
//
T_bpcfloat phase__in[NFREQS], T_hpcfloat **realj5t, //
T_hpcfloat
pt, T_hpcfloat freq[NFREQS],
//
int initx, int inity, int endx, int endy)
eise ii( m = 2 )
{ ^ = fopen("amp_iinc3y","w");
fpphase = fopen("phase_luie3y”,"w");
>
{
d i f f - endy - inity;
foi(i=inity;j<=inity+difpj++)
int i,j,m;
FILE
*fpphase;
T_hpcfloat phase_rad, phase_deg;
int i f f ;
{
amp[init![][j] = (l.O/amp_m[m]) *
sqrt(pow(real_pt[raitx]0],2.0) + //
pow (im agj)t[initx]0],2,0));
/*
% 9,8f \n",amp[iai£x][j]); */
fprintf(%," % 10.9f \n",amp[initx]OD;
/* phase[initx]0] = atan2(in3ag_pt[initx]jj],realjt[initx][j]) - U
phase_in[m];
*/
phase[initx][j] = atan2(imag_pt[initx][j],rcalj)t[initx][i]);
/* Save the phase values in range 0 to -360 */
phase_rad = phase[initx][j];
1) /* Along X */
{
for(m=0;m<NFREQS;m-H-)
{
if (m = 0 )
{ fp = fopen("anp_linelx'V’w");
^ ph ase = foperi(‘’phase__line!x'V'w");
}
else i f ( m ~ l )
{
= fopen(*'amp_lme2x","w");
fpphase = fopen("phase_line2x","w");
if(phase[imtx][j]>0.0)
{
phase_rad = phase[initx]0] - 2®Pi;
}
else i^m ==2)
{ = fopen("amp_line3x”,"w''};
^ p h ase = fopen("phase_line3x'V’w");
/* Converting phase to degrees. */
phase_deg = phase_rad*( 180/PI);
/* iprmd(^phase," % 9.8f W ,phase_deg}; */
fprintf^^phase," % 10 .9 f \n",phase_deg);
}
}
d ilf = endx - initx;
for(i=mitx; i<==iaitx+diff; i-H-)
fcto se® );
137
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
char dummy;
T_hpcfloat aux_f!oat;
fdose(fpphase);
} /» End o f FOR®/
}
}
fit = fopen("Ga","r");
Save the value o f Ez field along a line, wiiich can be either in
X or Y axis,
axis: i-> X , 2->Y
for(j=0;j<JEy++)
I
for(i=0;i<IE;i++)
{
fecanf(lp,"%f',&aux_float);
ga[i][j] = aux_flo 2 t;
Save_Ez_Iinex_wr{FILE ®fpEziinsxwr, T_hpcfioat **ez, int axis, //
int itiitx, int inity, int endx, int endy)
{ int i, j, difference;
T_hpcfloat value;
}
}
fclose(fp);
}$ 4c4c$4s4i4 : 4c4=4c 4:4=4:« 4c4c4i4 : 4c4: 4c4:4s4c4:4c4c4:4c $ 4c4i4>4c4c 4c4c4:4:4c4:4i4c
fscanf(fp,"%c",&dummy);
ii^axis— 1) /® A long X ® /
{
difference = endx - initx;
!St4:4ciSe4:s!!**4c4:4c!jc4c4t4c4s*«s4i*4t#!S4s4c4s4c4!45*4cS!4c*4:4=>!s*4c4s4<iSs4<*4c*s^!4;#**$*y
for(i=initx;i<=imtx+difFerence;i++)
{
Read_Gb(T_hpcf!oat gb[IE][JE])
value = ez[i][inity];
/ • 4>rintf(%," % 9.8f\n”,value); • /
/* %rintf(i)jEzlinexwr," % i5 .! 4 f ’,value); ®/
fprintf(fpEztinexwi," % 9.8f,value);
{
fprintfl^^Ezlmexwr,"ln");
FILE
in tij;
char diminiy;
T_hpcfloat aux_float;
if(a x is = 2 ) /* Along Y ®/
fp = fopen("Gb",'’r");
for(j=Oj<JE;j-H-)
}
{
for(i=0;i<IE;i++)
{
{
difference = endy - inity;
for(j=mity3 '<=inity+difference;j++)
{
fscanf{fp,"%f',&aux_float);
gb[i][j] “ aux_float;
value = ez[initx][j];
/ » fprmtfi;fi5,” %9.8f\n",value); */
/* iprintf(fi)Ezlinexwr,'' % 15.14f,value); */
iprintf{lpEzlinexwr,'' %9.8f',value);
}
}
fc lo se ® ;
}
iscanf(§>,"%c",&dutmny);
}
fprintfl[fpEz!inexwr,"'\n'');
}
}
ij:« « 4c
#4:d! jcti
^4c4t^4sj}:!jt*:S4i*4ciSe>!«#:*4c4c4s4t45 4s4s*4:S:4c4c4s4:>i(4c4s*4i##St4eSsi>4c*>is4!Jit4s4ccS454!S:!S4!
FuncEon Read_Parameters
This fimction reads the initiaiizing parameters from a file.
«#««i4c4c4c4s4:4c«4t4c««4c«c3'«»««Nc4c4c4c](84-.4s4c4c4s4>«4:4c4c4s«3!t4c4c4c4:««c«««c«4c»4cy
^#4:1$:
^
^ ^
^
^ >j!
void
Read_Paxameters(T_hpcfloat ’^dx,T_hpcfloat
*source_frequency,T_b 4 )cfloat *tO, //
T_hpcfioat *spread,int *npml,int *movie,int *movie_period, //
int "^tepSjint ’^axis,itit ’^initx,int *inity,int *endx,int *endy)
{ FILE
char variabie_name[80], movie_answer[ 15];
char file_name[80] = "inpiit.txt", dummy;
T_hpcfloat aux_float;
int index, aux_int;
read_param.h
void Read_Ga(T_hpcfioat ga[IE][JE]);
void Read_Gb(T_hpcfloat gb[IE][JE]);
void Read_Parameters(T_bpcfloat *dx, T_hpcfloat
*source_irequency, //
T_hpcSoat ’HO, T_hpcfioat ^spread, int ^ p m l, int ’^movie, //
int *moviej5eriod, int *nsteps, int *axis, int ’^initx, //
int
int *endx, int ^endy);
void Read_Properties(T_hpcfioat g[IB][JE], char filenameg[80]);
void Read_Properties_Fileaames(cbar filenamega[80], char
filenamegb[80j);
void Read__Steady_State(mt *last__period);
void Read_ID_fiicnaiT!e(char filenameID[SO]);
void Read_Prop_ID(int **proplD, char fiienameID[80]);
void Read_Only_Ezline(int *only_Ezlinexwr);
void Read_Start_Savmg^Ezlinex(int *start_saving_Ezlinex);
= fopen(f!!e_name,"r"))==NULL)
{
printfi[|''Can’t open file %s\a",file__name);
exit(0);
while(strcnp(vanabie_name,"START_PARAMETERS"))
fscan^fipmput,"%s",\^iab!e__Qame);
void
Read__Ga(T_hpcfloat ga[IE][JE])
fecanf(fpniput,”%s",variable__name);
frcanf(^input,"%s'',variabie_iiame);
fecanf(^input,"%f',dx):
{
HLE
int i,j;
/* DIM ENSIONS */
138
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
fscanf(fpmput,"%s",variable_nanie);
SOURCE
fecanf(^mput,"%s",v£riab]e_name);
exit(0);
}
V
fscanf(^input,"% f’,source_frequfincy);
fscan^45mpiit,''%s",variable_iiame);
fscanf(4^input,"%f,iO);
fscanf(^input,"%s",variabie_iiame);
fscan^^input,”% f’,spr©ad);
}
fscanf(^inpiit,"%s”,variable_naine);
fscanfl^^mput,"%s",variable_name);
fscaBi(fyinput,''%d'',npml);
/* PML «/
fscan^^input,'’%s",variable_name);
fscanf(^mput,"%s",variab!e_naine);
fscanfl|^^input,"%s",movie_answer);
/♦ M O V IE */
void
Read Properties{T_hpcfloat gPE][JEj, char fi!enaraeg[80])
{
F IL E *^ ;
in tij;
char dummy;
T Jipcfloat aux_float;
i6[lstrcmp(movie_answer/'yes'’))
if({% = fopen(Slenam eg,"r"))=NULL)
*movie = YES;
{
prmtfC'Can’t open file %s\n",filenameg);
exit(0);
else if(lstrcnip(mo^de_ajiswer,'*i2o”))
}
^movie = NO;
for(j=Oj<tEj-H-)
eise
printf("Undefmed answer to movie W );
exit(0);
{
for(i=0;i<IE;i-H-)
{
fscanf(fp,''%f',&aux_float);
g[i][j] = aux_float;
fscanf(fpii^ut,"%s",variable_name);
fscanf(f^input5”%d",movie_period);
}
}
fclose(^ );
}
fscanfl;^,"%c'',&damniy);
fscanf(^input,"% s",wiable_name);
fecanf(^input,"%s’',variabie_name);
fscanfl^^input,''%d",nsteps);
/* ITERATIONS */
fscaj3f[^input,"%s";yariable_name);
fscani(^mput,"%s",,variable_name);
fscan^^mput,"%d”,axis);
fecaaf(^input,"%s";,variable_name);
&caef(^mpiit,"%d’ ,initx);
fscanf(^mput,"%s".',variable_name);
fscaid(^m put,”%d’',inity);
fscanf(Q)raput,''%s"',variable__narae);
fscanf(^input,”%d’',endx);
fecanf(^mput,"%s"'»variable_Qaine);
fscanf(^mput,"%d’',endy);
/ • SAVE LINE • /
Function Read_Properties_Filenames
This fimction reads the names o f the files containing properties Ga
and Gb.
void
Read_Properties_Filenames(cliar fl!enamega[80], char
filenamegb[80])
{ FILE "^^input;
char variable__name[80];
char Sle__name[80] = ''input.txt", dummy;
if((^m put = fopen(fiie__name/'r"))— NULL)
{
fclose(fpinput);
printf("Can*t open file %s\n",file_name);
exit(0);
printf(”ddx is % 8.7f W’,*dx);
print^"Freqiiency is %2. I f GHz \n'',’^source_fequeiicy/le9);
printf("tO is % f \n",*tO);
printf("spread is % f \n",*spiead);
print^"npml is %d \n",^npmi);
print^''movie is %d \n'',*movie);
print£("movie_period is %d \n",*moviej>eriod);
printi(”nsteps is %d \n”,*r^teps);
printS[^''axis is %d \n",*axis);
printfl;"initx is %d \n",^initx);
printf("imty is %d \n",*mity);
printf("emix is %d V ,*en d x);
printf("endy is %d \n‘',®endy);
while(strcmp(variable_name,"PROPERTIES_FILES.-"))
fscanf(fpinput,"%s",variable_name);
lscanf(^mput>"%s",variable__name);
6caaf(^input,"% s”,filenamega);
fscanf(^input,"%s",variable_name);
fscani(^input,''%s”,filenamegb);
fclose(^input);
printf|‘'Ga file is %s \n",filenamega);
printfl^^Gb file is %s \n'',filenamegb);
/* GA: */
/* GB */
printf("Are parameters correct? (y /n ):");
scanf("%c",&dummy);
if(dunim y='n’)
139
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Function Read_Steady_State
for(i=0;j<JEy++)
(
for(i=0;i<IE;i-4-4-)
{
fscan^;4>,"%d",&aux_int);
void
Read__Steady State(int '^iastjjeriod)
(H L E ^ fp ;
char variab!e_name[8GJ;
char fiie_name[80] = "iBput,txt"j dummy;
propID[i][j] = ai!x_int;
}
fscanf(^,"%c",&;dummy);
}
fclose(fb);
}
il((lp = fop^n(fi!e_Bams,"r”))==NULL)
{
printf("Can't open file %s\n",f!le_name);
exit(0);
Function Read_Oniy_Ezlme
while(strcmp{variable_tmme,"STEAD Y_STATE.-"))
{
}
void
Read Oniy_Ezline(mt *only_Ezlinexvv'r)
i,"%s",variable__nan:K);
{
FILE^^;
fscanf(^,"%s",variable_name);
/* Lastjperiod: */
fscanf[^,"%d",last_period);
print^"Last period is %d ^B'’,*lastj3criod);
fciose(^ );
char variable_narne[80], ezliiiexwr_answer[ 15];
char fiie_name[80] = "input.txt", dummy;
if((^ = fopcn(fiie_name/'r"))==NULL)
{
}
prmtfC'Can’t open file %s\n",file_name);
exit(0);
}
Function Read_ID_Fiiename
This fimction reads the name o f the file containing properties ID.
whi!e(strcn^(variable_name,''OPTION__OUTPUT.-"))
{
fscanfi^^,"%s",variable_name);
}
Read_ID_jFilename(char filename!D[80])
{ FILE *^input;
char variable_name[80];
char file_name[80] = "input.txt", dmnmy;
fscanf(fp,”%s'',variable_name);
fscanf(^,"%s",ezlmexwr_answer);
/* only__E2 linexwr:
if((^m put = fopen(fils name,"r"))==NULL)
(
if(!strcmp(ezlinexwr_ai^wer,"yes"))
printf("Can't open file %s\n'’,file_name);
exit(O);
*oniy_Eziinexwr = YES;
}
else ifi^!strcn^(ezlinexwr_answer,"no"))
while(strcmp(variabk_narae/'PROPERTIES_FILE_ID.-"))
I
*only_Eziinexwr = NO;
}
fscanf(^!nput,"%s",variable_name);
fscanf(^inpiit,"%s",mriable_name);
fscanf(^input,"%s",fiienameID);
fclose(^input);
prmtf["Prop__ID file is %s \n",filenameID);
else
{ printfC'Undefined answer to only_Ezlinexwr \n");
exit(0):
PropJD : */
}
prmtfC’Only Ezlinexwr is %d \n",*only_Ezliiiexwr);
fclose(%);
return;
Function Read_Start_SaYing_Ez!inex
void
Read_Prop_ID(int **proplD, char fi!enameID[80])
void
Read_Start_Saving_Ezlmex(int *start_saving_Ezlinex)
{ FILE
char variab!e_name[80];
char file_name[80] = '’input.txt", dummy;
{
FILE»^;
int ij;
char dummy;
int aux_mt;
il( ( ^ ~ fopen(file_nam e,"r"))=N ljLL)
{
= fopen(fdenanieID,"r’'))==NULL)
{
prmtf("Can't open file %sNn"/ile_name);
exit(0);
prmtf("Can’t open file %s\n",fiienameID);
exit(0);
140
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
while(strcit 5 s(vanabie_name,"START_SAVING_EZLINE.-"))
T_hpcfloat ga[IE][JE], T_hpcfioat gb[IE][JE], TJipcfloat
iz[IE][JEJ. //
T__hpcfloat freq[NFREQS3, T_hpcfloat arg[I^JFREQS], T_hpcfloat
n;aTjn[NFREQS], //
T__hpclIoat imag_m[NFREQS], T_hpcfloat amp__in[NFREQSj, //
T__hpcfloat phase_in[NFREQS], T_hpcfloat
reai_pl[NFREQS3[IE3[JE3, //
T_hpcfloat im agjit|NFREQS3[IE][JE], T___hpcfloat arapflEjfJE],
T_hpcfioat phase[IE3[JEJ)
{
fscan^fp,"%s",vsrmble_name);
fscanf(fp,"%s",variabie__nanie);
/* start_saving_Bzlmex: */
fscan^^,"%d'’,stait_savmg_Ez!mex);
printf("Start Saving Ezlinex Stq? is %d \n",*start_saving_Ezlinex);
fciose(^ );
return;
{
int y ,m ;
}
for(j=Oa<JEy++)
{
ez_mc{j] = 0.0;
h x jn cfj] = 0.0;
initiaiize.h
fo^i=0;i<IE;!++)
{
void Initialize ArraysCT^hpcfloat ez_inc[JE] 5T_hpcfloat hx_inc[JE],
//
T_hpcfloat ez[IE][IE],T_hpcfloat dz[iE][JE],Tj 0 pcfloat
hx[IE][JE], //
T_hpcfloat hy[IE][JE],T_hpcflost iiiK;[IE][JE|,T_hpcfloat
ihy[IE][JB],//
T__hpcfloat ga[IE][JE],T_hpcfloat gb[IE][JE3,T_hpcfloat
iz[lE][JE], //
TJypcflo& t fi«q[NFREQS],T_hpcfloat arg[NFREQS],
//
T_hpcfloat real__in[NFREQS], T_hpcfloat imag_m[NFREQS], //
T_hpcfIoat amp_in[NFREQS], T_hpcfloat phase_m[NFREQS],
ez[i][j] = 0.0
dz[i]03 = 0.
ltx[i][i] = 0.
hy[i][13 = 0.
ihx[i]0J = O
ihy[i30] = 0
g a [i]D ]= l.
gb[i][j] = 0.
iz[i]{j] = 0.0
//
T_hpcfloat reai_j)t[NFREQS][IE3[IE3, T__hpcfloat
ima^_pt[NFREQS][IE][JE], //
T_hpcfloat axnp[IE][JE], T Jipcfloat phase[IE][JE3);
for{m=0;m<NFREQS;m-H-)
(
void Initialize_Arrays2(T_hpcfloat ez__inc[JE],T_hpcfloat
hxJnc[JE],
//
T_hpcf!oat ez[IE][JE],T_lipcfloat dz[IE][JE],T_hpcfloat
hx[IE][JE3, //
T_hpcfloat hy[IE][JE],T__hpcfloat ihx[IE][IE],T_!:^fIoat
ihy[lE3[JE],//
T_hpcfloat ga[IE3[JEJ,TJipcfloat gb[IE][JE], T_hpcfloat
izpE3[JEj);
void Initialize_AiTays3(T___bpcfloat ez_inc[JE3, TJipcfloat
hx_inc[JE], //
T Jipcfloat ^^ez, T_hpcfloat
T__I^float **hx,
T_i^cfloat ^*hy, T_b|x:float **ihx, T_hpcf!oat
TJipcfloat ^*iz);
/* gb = (cond * dt) / epsz */
/* iz - iz + gb*Ez
freq[m3 = 0.0;
aig[m] = 0.0;
realjn[m ] = 0.0;
im agjn [m ] = 0.0;
amp__m[m] = 0.0;
phase_^mj =0.0;
for(j=Oy<rEy++)
{
for(i=^;i<IE;i-H-)
{
real_pt[m][i][j] = 0.0;
«na£J>t[tn][i]Li]=0.0;
//
//
an^[i][j] = 0.0;
phase[i][j] = 0.0;
)
void initiaIi2 e_A rraysJD (int **prop_ID, //
T_hpcfloat constGa[TOMBER_OF__TTSSUES3, T_hpcfioat
cofistGb[NUMBER__OF_TISSUES3);
void initialize_Arrays_Fourier(T_hpcfloat **i'eal_j>t, T_hpcfloat
**imag_pt, //
T_i^cf!oat ^*amp5 T_hpcf!oat ^ ^ h ase,
//
T _l:^ float real_m[NFREQS], T_hpcfloat
imag__iii[NFREQS3s //
T_hpcSoat amp__in[NFREQS], T__hpcfloat
phase_inp'^FREQS], //
T_hpcfloat freq[NFREQS], T_hpcfloat arg[NFREQS]);
void
Initialize_Arrays2(T_hpcfloat ez_inc[JE], T_hpcfioat hx_inc[JE], //
T_bpcfloat ez[IE][JE], T_hpcfloat dz[IE][JB], T_bpcfloat
hx[IE][JE], //
T_hpcfloat hy[IE][JE], T_hpcfloat ihx[IE][JE], T_hpcfloat
ihy[!E][JE], //
T_hpcfloat ga[IEj[JE], T_hpcfloat gfa[IE][JE], T_hpofloat
izpE][JE])
{
Iratiali2 e_Arrsys(T_hpcfloat ez_inc[JE], T_hpcfloat hxJnc[JE], //
T__hpcfl[oat ez[!E][JE3,T_bpcfioat dz[IE][IE], T_hpcf!oat
hx[IE3[JE],
T__hpcfloat hy[IE][JE], T_hpcf!oat ilix[IE3[JE|, T_lipcfloat
ihy[IE3[JE3, //
in tij;
n
for(j=Oa<JEy++)
{
ez_inc[j] = 0.0;
141
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
hx_itic[jj = 0.0;
void Initialize_Arrays__Fouri€r(T_hpcfioat ^*real_pt, T_hpcfioat
for(j=0;i<IE;i++)
**iinagj5t, //
T_hpcfIoat *^amp, T_hpcfloat ’^‘^phase,
//
T__hpcf3oat r e a H ti[N F ^ Q S ], T_hpcfloat
im agJn[NFREQS], //
T_hpcfl.oat amp_in[NFREQS], T_lipcf!Gat
pb.ase_in.[NFREQS], //
T_}pcfloat Ireq^NFREQS], T Jipcfloat arg[NFREQS])
ez[i][j] = 0.0
dz[i]D] = 0.0
MOD]=0-0
hy[i]0] = O.O
ihx[i]0] = O-'
ihy[i]0] = O,i
g a [i]0 1 = l.O
gb[i][j] = 0,0
iz[i]0] = O.O;
{
int ij,m;
gb = (cond * dt) / epsz
/ • iz = iz + gb*Ez
*!
for(fTF=0;ni<NFREQS;raH-t-)
{
&eq[m] = 0.0;
arg[m] = 0,0;
reai_m[m| = 0.0;
mmg_in[m] = 0.0;
amp_in[m] = 0.0;
ph^e_in[m ] = 0.0;
void
Imtialize_AiTays3(T_hpcfloat ez_mc[JE],T_hpcfloat ta_mc[JE], //
T 3 p o flo a t **ez, T_hpcfloat **dz, T_hpcfloat **hx, //
T_hpcfloat **hy, T__bpcf1ost **ibx, T^hpcfloat ®*siiy, //
T_hpcfloat **iz)
for(j=0;J<JE;j++)
{
for(i=0;i<IE;i++)
{
real_pt[i]0] = 0.0;
Di'ag_pt[i]0]=0.0;
{
in tij;
for(j=O yO Ej++)
amp[i]|J3 ~ 0.0;
phase[i][jj = 0.0;
{
ez_inc[j] = 0.0;
iix_inc|j] - 0.0;
}
for(i=0;i<IE;i++)
{
ez[i][j] = 0.0;
dz[i][j] = 0.0;
hx[i][j] = 0.0;
hy[i]0] = O.O;
ihx[i][j] = 0.0;
ihy[i][j] = 0,0;
iz[i]0] = 0.0;
an:^Ez_ssi^c.h
void Save_ArapEz(T_hpcfloat ^*amp_ez);
T_hpcfloat Mag(T_hpcfloat valne);
void Save_An]pEz_Lme(mt axis, T_hpcfloat ’^^arap_ez, int initx, int
in ity ,//
int endx, int endy);
void Reset__AnfiEz(T__hpcf!oat **arap_ez);
void Update_AinpEz(T_hpcf!oat
T_l^cfloat '*=*an^__ez);
/* iz = iz + gb*Ez “/
}
}
void
Initia!ize_Arrays_ID(int ®*propID, T_hpofloat
constGa[NUM BER_OF_TISSUES], //
T_lipefioat constGb[NUMBER OF_TiSSUES])
Write the amplitude o f Ez field out to a file
{,
{
"AnpEz"
void
Save_AmpEz(T_l3pcfloat *^anip__ez)
,
int ij;
char filenajtie[20];
HLE*^;
foi(i=O j< JE j-H -)
{
ifl^(55=fopen("AtnpEz","w"))=NULL)
{ p rin tftjiT h e file AmpEz could not be opened");
exit(0);
for(i=0;i<IE;i++)
{
propID[i]0] = FREE_SPACE;
}
for(j=OjOEj++)
for(i=0;i<NUMBER_OF_TISSUES;i+y)
{
{
fof(i=0;i<IE;i++)
{
ijsrintfl'lp,” % !0.9f',ainp_ez[i][j]);
}
^rintf(ip," \n");
}
cons£Ga[i] = 1.0;
constGbp] = 0.0;
142
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
fprintf(ii>," % 10.9f W .vatoe);
fclose(lp);
}
}
fclose(fp);
Update the atnplitade o f E z fieid.
}
R eset the array containing the amplitude o f E z fieid
Update_AmpEz(T__hpcfloat ’^“^ez, T_bpcfloat "®*anip_ez)
{
void
Reset_AmpEz(T_hpcfloat **atiip_ez)
in tij;
T_hpofloat roag_Ez;
{
for{j=Oj<JEy++)
in tij;
{
for(i=0;i<iE;i++)
{
for{j=Oj<JEj++)
I
ioi(i=fl;i<IE;iH-+)
mag_Ez = M ag(ez[i]0]);
if (mag_Ez > am p_ez[i]0])
amp_ez[i][j] = m a g jiz ;
(
amp_ez[i][j] = 0,0;
}
>
}
}
I
T__hpcfloat
Mag(T_hpcfloat value)
nKm_n:]gpit.h
(
T_hpcfloat *’’'Res_Mem_T_hpcfloat(mt max__x,ini max__y);
int '*^Res_Mera_int(int max_x,mt max_y);
void Free_Mem_T_hpcfloa^T_fapcfloat **array, int max_x, int
raax_y);
void Free_Mem_int(int •'•'array, int inax_x, int max_y);
T_hpcfloat iraa^itude;
if(value < 0.0)
m ^ t u d e = value * (-1);
else
nrnguitude = value;
retum(magnitude);
7
}
T_hpcfloat
**Res_Mem_T__hpcfioat(int raax_x, int rmx_y)
Save the amplitiKle o f the Ez field along a line, which can be either
in
X or Y axis,
axis: 1->X, 2->Y
{
inti;
T_hpcfloat *’*‘array;
if?t(array = (T__hpcf!oat **)malloc(sizeof(T_hpcfIoat *)*nmx_s)))
void
Save_AmpEz Lme(int axis, T_hpcfloat **ampEz, int initx, int inity,
{
//
printS;"Cannot allocate T_hpcfloat raemory\n");
exit(l);
int endx, int endy)
(
}
int i, j, difference;
FILE
T_hpcfioat value;
ifl[axis==l)
for(i=0;i<inax_x;i-H-)
array[ij = (T__hpcfioat *)calloc(im x_y, sizeod^T_hpcfloat));
retum(array);
Along X * /
}
{
^ = fopen(”ampEz_linex'V’W');
difference = endx - initx;
for(i=initx;i<=initx+difference;i++)
{
int
**Res_Mem_int(int max_x, int max_y)
value = anpEz[i][inityj;
^rint^%," %10.9f\n",value);
{
}
fclose(^ );
}
int i, **array;
ifl;!(array = (int **)maUoc(sb:eo^mt ‘^)*imx_x)))
{
if(a x is = 2 ) /* Along Y */
printfC'Cannot allocate INT meraory\n");
exit(l);
{
fp = fopen("an^E 2 _liney’V'w");
for(i=0;!<max_x;i-H')
array[i] = (int *)caJloc(max__y, sizeof(int));
difference = endy - inity;
for(j=inityj<=inity+difference'J++)
{
retum(array);
value = ampEz[initx]jjj;
143
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
)
i^(^ fopen(" ave_am p E z" ," w " ))=N U L L )
( prinifl^"\nFile A ve_A np E z conid not be opened”);
exit(0);
void
Free_MeRi_T_hpcfloat(T_hpciloat **array, int max_x, int m ax_j)
for(J=0;j<JEy-H-)
int i;
{
for(i=0;i<IE;i-H-)
{
}
for(i=0;i<inax_x;i++)
free(airay[i]);
% 15.14f',ave_am pez[i][j]);
tjprintfl^fp," \n");
free(array);
}
)
fclose(^ );
}
Update the averaged amplitude o f E z field.
pn = Current Ez value.
p n l = Previous Ez value.
pn2 = Two time steps back Ez value.
void
Free_Mem_int(int ’’'^array, int max_X5 int mas_y)
{
inti;
void
Update_Ave_AmpEz(T_hpcfloat *'^ez, T_hpcfloat **ave_ampez, //
T_hpcfloat **one_stepJ)ack_ez,
T_hpcfloat **two_step_back_ez,
//
int **peaks_count8r)
{ in tij;
T Jipcfloat pn, p n l, pn2;
T_hpcfloat average, prod;
for(i=0;i<max_x;i-H-)
free(array[i]);
free(array);
)
arapEz_sshpc.h
forO=OJ<JEJ++)
{
void Copy_fArrays(T_hpcfloat **source, T_hpcfloat ^‘^target);
void Save_Ave_An^Ez(T_hpcfloat **ave_ampez);
void Save_Ave__A2npEz_Line(int axis, T_hpcfloat **ave_ampe 2 , int
initx, //
int inity, int endx, int endy);
void Reset_fArray(T_bpcfloat ^*ave_ampez);
void Reset_iArray(int **peaks_counter);
void Update_Ave_An:roEz(T_hpcfloat *^ez, T_l:^ciloat
**ave_arapez,
//
T_hpcfloat **one_step_back_ez,
//
T_hpcfloat **two_step_back_ez,
//
int **peaks_coimter);
void Read_Average_Steady_State(int *start_averaging_step);
TJipcfloat Absolute(T_hpcfioat value);
for(i=0;i<IE;i-^+)
{
pn = ez[i][j];
pnl = one_step_back_e 2 [i][j];
pn2 = two_step_back_ez[i][j];
/* p nl is a positive peak p n l is a negative peak */
ii( ((p n l> p n )& & (p n l^ n 2 )) jj ((pni<pn)& & (pnl <pn2)) )
(
average = ave_ampez[i]{j];
prod == average * peaks_counter[i][j];
-H-peaks_counter[i][j];
ave_ampez[i][j] = ( prod + A b so lu te{p n l)) /
peaks_coimter[i] [j ];
void
Copy_fl\rrays(T_hpcfloat ^^source, T__hpcfloat **target)
{ int ij;
ii((i= IC )& & (j= = 6 1 0 ))
{
printf(''\nPeak
");
printf("PeaJc value = % 9.Sf ",pnl);
printf("Peaks counter = %dW,peaks_counter[i][j]);
for(j==Ou<JE;j-i-+)
{
}
for(i=0;i<IE;i++)
*/
}
target[i][j] = source[i][j];
}
}
T_hpcSoat
Absohite(TJipcfloat value)
{ T_bpcfloat absolute;
void
Save_Ave_AmpEz(T_hpcfloat ®®ave__ainpez)
{ int U;
char fi!eimme[20];
FILE *fy ;
if(value < 0.0)
absolute - value ’ (-1);
else
144
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
absolute = value;
retum(absoiute);
}
Function Read_Steady_State
Save the amplitude o f the Ave_AmpEz field along a line, which can
be either in
X or Y axis,
axis: 1->X, 2->Y
void
Read_Average_Steady_State(int "^step)
{ FILE '^fp;
char variable_name[80];
char file_narae[20] = "input.txt", dummy;
void
Save_Ave__AmpEz_Lme(int axis, TJipcfloat *®ave_anipez, int initx,
int inity, //
int endx, int endy)
{ int i, j, difference;
FILE *fp;
T_hpcfioat value;
if((ip = fopen(file__name,"r")}==NULL)
(
printf("Can't open Eie %s\n",file_name);
exitfO);
}
if ( a x is = i ) I* Along X
{
while(strcn^(variable_name,"AVERAGE__STEADY_STATE.-"))
{
^ = fopen("ave_ampEzJmex","w");
difference = endx - initx;
for(i=mitx;!<=initx+difference;i-K-)
ffscanf(fp,'’%s",variable_name);
}
/* start_averaging^step V
fscanf(fp,"Yos",variable__name);
fk;anf(^,"%d’',step);
printf("Start averaging step is %d '®.'',*step);
(
value = ave__anpe2 [i][inity];
/* fprintf(^,’' %9.8f\n",value);
^rintf(^," %15-14f\n",value);
fclose(^ );
}
fclose(^ );
}
}
ifl^axis=2) /* Along Y ®/
{
^ == fopen('’ave_ampE 2 _liney'V'w");
difference = endy - inity;
for(j=inity;j<=Tmty-5-difference;j-H-)
{
value = ave_ampez[initx]0];
“/o9.8f\n",value); ’*'/
^rintf(fy," %15.14f\n",value);
1
fclose(^ );
}
}
R eset a float array
void
Reset_fArray(T_hpcfIoat **array)
{ int iJ;
for(j=Oy<JE;j++)
{
for(i=0;i<IE;i++)
{
array[i][j] = 0.0;
}
}
Reset an int array
void
Reset_iAiray(int '**array)
{ int g ;
for(j=Oy<JEa-H-)
{
for(i=0;i<!E;i-H-)
I
array[i][j] = 0;
145
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
REFERENCES
[1] American Lung Association of Hudson Valiey
http://www.alahv.org/bookfiles4/iung_cancer.Mml
[2] P. Kumar, M. Clark, “Clinical medicine: a textbook for medical students and
doctors”, 4* ed. Edinburgh: W.B. Saunders, 1998.
[3] L. E. Larsen and J. H. Jacobi, “Medical applications of microwave imaging”, New
York: IEEE, 1986.
[4] E.C. Fear, S.C. Hagness, P.M. Meaney, M. Okoniewski and M.A. Stuchly, “Nearfield imaging for breast tumor detection”, IEEE Microwave Magazine, vol. 3, March
2002, pp. 48-56.
[5] P. M. Meaney, K. D. Paulsen, A. Hartov, and R. K. Crane, “An active microwave
imaging system for reconstruction of 2-D electrical property distributions”, IEEE
Transactions on biomedical engineering, vol. 42, No. 10, October 1995, pp. 1017-1026.
[6] Fear E.C., Meaney P.M., and Stuchly M.A., “Microwaves for breast cancer
detection?, IEEE Potentials, vol. 22, Issue 1, Feb. - March 2003, pp. 12-18.
[7] S. C. Hagness, A. Taflove, and J. E. Bridges, “FDTD Analysis of pulsed microwave
confocal system for breast cancer detection”, Proceedings -
19* International
Conference - lEEE/EMBS, Oct. 30 - Nov. 2, 1997 Chicago, IL. USA
[8] S. C. Hagness, A. Taflove, and J. E. Bridges, “Three-dimensional FDTD analysis of
a pulsed microwave confocal system for breast cancer detection: Design of an antennaarray element”, IEEE Transactions of antennas and propagation, vol. 47, No. 5, May
1999, pp. 783-791.
146
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[9] X. Li, and S. C. Hagness, “A confocal microwave imaging algorithm for breast
cancer detection”, IEEE Microwave and wireless components letters, vol. 11, No. 3,
March 2001, pp. 130-132.
[10] E. C. Fear, X. Li, S. C. Hagness, and M. A. Stuchly, “Confocal microwave imaging
for breast cancer detection: Localization of tumors in three dimensions”, IEEE
Transactions on biomedical engineering, Vol. 49, No. 8, August 2002, pp. 812-822.
[11] X. Li, S. C. Hagness, B. D. Van Veen, and D. van der Weide, “Experimental
investigation of microwave imaging via space-time beamforming for breast cancer
detection”, 2003 IEEE MTT-S Memational Microwave Symposium Digest.
[12] M. Kubo, Y. Kawata, N. Niki, K. Eguchi, H. Ohmatsu, R. Kakinuma, M. Kaneko,
M. Kusumoto, N. Moriyama, K. Mori, and H. Nishiyama, “Automatic extraction of
pulmonary fissures from multidetector-row CT images”, International Conference on
Image processing, 2001, vol. 3, October 2001, pp. 1091-1094.
[13] C. B. Caldwell, K. Mah, Y. C. Ung, C. E. Danjoux, and J. M. Balogh, “FDGPET/CT Integration: Impact on tumor localization an dose volume histograms in
radiation therapy”, IEEE Procedings of the 22“^*annual EMBS international conference,
July 23-28, Chicago, IL.
[14] M. N. O. Sadiku, and A. F. Peterson, “A comparison of numerical methods for
computing electromagnetics fields”, Southeastcon ’90, Proceedings, IEEE , Vol. 1, 1-4
April 1990, pp. 42-47.
[15] K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell’s equations in isotropic media”, IEEE Transactions on antennas and
propagation. Vol. AP-14, May 1966, pp. 302-307.
[16] A. Taflove, and M. E. Brodwin, “Numerical solution of steady-state
electromagnetic scattering problems using the time-dependent maxwell’s equations”,
147
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
IEEE Transactions on microwave theory and techniques^ Vol. MIT-23, No. 8, August
1975, p p . 6 2 3 -6 3 0 .
[17] A. Taflove, M. J. Piket-May, and S. C. Hagness, “FDTD - How complex a
problem can we solve”, Antennas and propagation society international symposium,
2000. IEEE, Vol. 3, 16-21 July 2000, pp. 1641.
[18] Shlager, K. L., and J. B. Schneider, “A selective survey of thefinite-difference
time-domain literature”. Antennas and propagation magazine, IEEE, Vol. 37, Issue 4,
August 1995, pp. 39-57.
[19] A. Taflove and M. E. Brodwin, “Computation of the electromagnetic fields and
induced temperatures within a mode! of the microwave-irradiated human eye,” IEEE
Trans. Microwave Theory Tech., vol. 23, pp. 888-896,1975.
[20]
Web
page
of United
Kingdom
National
CancerInformation
service
[http://www.cancerbacup.org.uk/]
[21] T. Okumura, T. Miwa, J. Kako, S. Yamamoto, M. Matsumoto, Y. Tateno, T.
linuma, and T. Matsumoto, “Automatic detection of lung cancers in chest CT images by
variable N-Quoit filter”, Proceedings. Fourteenth International Conference on Pattern
Recognition, vol. 2, August 16-20 1998, pp. 1671-1673.
[22] Web page o f CancerCare and the Oncology Nursing Society.
[http ://www.lungcancer.org/]
[23] B. Jacob, “Handbook of medical imaging”, Bellingham, Wash., 2000.
[24] K. Kanazawa, M. Kubo, N. Niki, H. Satoh, and H. Ohmatsu, “Computer aided
diagnosis system o f lung cancer based on helical CT images”, IEEE Transactions on
biomedical engineering, vol. 42, No. 10, October 1995, pp. 1017-1026.
[25] H. K. Son, M. Yun, T. J. Jeon, D. O. Kim, H. J. Jung, J. D. Lee, H. S. Yoo, and H.
J. Kim, “ROC analysis of ordered subset expectation maximization and filtered back
148
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
projection teciinique for FDG-PET in lung cancer”, IEEE Transactions on nuclear
science, vol. 50, No. 1, Februaiy 2003, pp. 37-41.
[26] M. A. Brown, and R. C. Semelka, “MRI Basic principles and applications”, 2"^^
Ed., New York John Wiley and Sons, 1999.
[27] Kunz, K. S., Luebbers, R. J. “The Finite Difference Time Domain Method for
Electromagnetics, CRC Press, 1993.
[28] D. M. Sullivan, “Electromagnetic simulation using the FDTD method”, IEEE Press
Series on RF and Microwave Technology, 2000.
[29] A.F. Peterson, S.L. Ray, and R. Mittra, “Computational Methods for
Electromagnetics”, IEEE Press, 1998.
[30] A. Taflove, “Computational Electrodynamics, The Finite-Difference Time-Domain
Method”, Artech House, Inc., 1997.
[31] J. P. Berenguer, “A perfectly matched layers for the FDTD solution of wavestracture interaction problems”, IEEE Trans. Antennas and Propagation, vol. 44, No. 1,
Jan. 1996, pp. 110-117.
[32] D. M. Sullivan, “A simplified PML for use with FDTD method”, IEEE Microwave
and guided wave letters, vol. 6, no. 2, February 1996, pp.97-99.
[33] J. Represa, C. Pereira, M. Panizo,” A simple demonstration of numerical dispersion
under FDTD”, IEEE Transactions on Education, Vol. 40, No. 1, February 1997, pp. 98102.
[34] W. H. Hayt, “Engineering electromagnetics”, New York McGraw Hill, 1989, 5*
Ed., pp. 359.
[35] R. F. Harrington, “Time-harmonic electromagnetic fields”, New York McGraw
Hill, 1961.
149
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[36] C.A. Balanis, “Advanced Engineering Electromagnetics”, 2”‘^ Ed., John Wiley &
Sons, Inc., 1989.
[37] K. Umashankar, A. Taflove, “Computational Electromagnetics”, Boston: Artech
House, 1993.
[38] S, Caorsi, E. Bermani and A. Massa, “A fmite-element procedure based on a
boundary-value approach for th e ' evaluation o f the electromagnetic exposure in
biological phantoms”, IEEE Transactions on Microwave Theory and Techniques, Vol.
50, No. 10, October 2002.
[39]
Stanford
University
Medical
Media
&
Information
Technologies,
http ://summit. Stanford. edu/RESEA R C IW isibleHuman/THORAX/A_vm 1451 .jpg
[40] D. Panescu, J.G. Webster, R.A. Stratbucker, “Modeling current density
distributions during transcutaneous cardiac pacing”, IEEE Transactions on Biomedical
Engineering, Vol. 41 Issue 6 , June 1994, pp. 549 -555
[41] A.M. Morega, D. Mocanu, M. Morega, A. Stefan, “A dynamic model for cardiac
defibrillation”. Engineering in Medicine and Biology Society,
1996. Bridging
Disciplines for Biomedicine. Proceedings of the 18th Annual International Conference
of the IEEE, Vol. 5, 31 O ct-3 Nov. 1996, pp 2291 -2292.
[42] Institute for Applied Physics web page, http://niremf.ifac.cnr.it/tissprop/
[43] C. Gabriel, S. Gabriel and E. Corthout, “The dielectric properties of biological
tissues: I. Literature survey”, Phys. Med. Biol., Vol. 41, Nov. 1996, pp.2231-2249.
[44] S. Gabriel, R.W. Lau, and C. Gabriel, “The dielectric properties of biological
tissues: II. Measurements in the frequency range 10 Hz to 20 GHz”, Phys. Med. Biol.,
Vol. 41, Nov. 1996, pp.2251-2269.
150
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
[45] S. Gabriel, R.W. Lau, and C. Gabriel, “The dielectric properties of biological
tissues: IH. Parametric models for the dielectric spectrum of tissues”, Phys. Med. BioL,
Vol. 41, Nov. 1996, pp.2271-2293.
[46] “Harper Collins Illustrated Medical Dictionary” , 4* ed., Harper Resource Harper
Collins Publishers, 2001.
[47] S.C. Hagness, A. Taflove, and J. E. Bridges,“FDTD Modeling of a coherentaddition anteima array for early-stage detection of breast cancer”, Antennas and
Propagation Society International Symposium, 1998. IEEE , Vol. 2 , 21-26 June 1998
pp. 1220 -1223.
[48] S.C. DeMarco, G. Lazzi, Wentai Liu; J.D. Weiland; M.S. Humayun; “Computed
SAR and thermal elevation in a 0.25-mm 2-D model of the human eye and head in
response to an implanted retinal stimulator - part I: models and methods”, Antennas and
Propagation, IEEE Transactions on,Vol. 51 , Issue: 9 , Sep. 2003, pp. 2274 - 2285.
[49] S.C. DeMarco, G. Lazzi, Wentai Liu; J.D. Weiland; M.S. Humayun; “Computed
SAR and thermal elevation in a 0.25-mm 2-D model of the human eye and head in
response to an implanted retinal stimulator - part II: Results”, Antennas and
Propagation, IEEE Transactions on,Vol. 51 , Issue: 9 , Sep. 2003, pp. 2286 - 2295.
151
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
BIOGRAPHICAL INFORMATION
Luis M. Camacho-Velazquez was bom in Poza Rica, Mexico.
He eamed a
bachelor degree at Universidad Autonoma de Nuevo Leon, Mexico in 1988. On the
same university he eamed a Master of Science degree on 1991. From October 1991 to
March 1993 he conducted a research stay at Technische Hochschule Darmstadt,
Germany.
Since 1989 he has worked as a professor at Universidad Autonoma de Nuevo
Leon.
He is currently working toward a PhD degree at the Electrical Engineering
Department, University of Texas at Arlington.
During his stay as student at UTA, he has been appointed as Graduate Teaching
Assistant at the Electrical Engineering Department.
His main research interests are numerical techniques for electromagnetics, and
microwave imaging.
152
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Документ
Категория
Без категории
Просмотров
0
Размер файла
8 029 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа