# Contributions to cost reduction and sensitivity improvement of microwave breast cancer detection

код для вставкиСкачатьCONTRIBUTIONS TO COST REDUCTION AND SENSITIVITY IMPROVEMENT OF MICROWAVE BREAST CANCER DETECTION by Min Zhao A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) at the UNIVERSITY OF WISCONSIN-MADISON 2009 UMI Number: 3400124 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. UMT Dissertation Publishing UMI 3400124 Copyright 2010 by ProQuest LLC. All rights reserved. This edition of the work is protected against unauthorized copying under Title 17, United States Code. ProQuest LLC 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, Ml 48106-1346 o ID W O 3" O u Q. C o 3 3 O O o u ? c 3 ai **• (A IQ * 3 a IB o o 3 5' a (A (0 * * * * * * ** * * * * ft » » t ** * » * ** s J 3 CO IQ 3 01 < — o •a > Committee's Page. This page is not to be hand-written except for the signatures (Q 3 5" • • * 01 3 5* CQ C (A > c K) O O •O c </> 3 w. > *< c o O (D 3 o- <D o fl> a a. CD Q. at 2L O SL m x (0 3 O o I <v -< a u B» 2 o Q O D M a- 3 (D 3 CD CD >< en 3 CD •g3" I~ OT O - * CD If -» 8" f D " o • 9. O -* a. 3 5 0) 3- CD CD O s 3" CD 3" o o o CO 0) c 0) a. Q 5 o C O ~ ccr C 3 3 3: <' QCD. in < —^ o o CO O C D Q D CO CD CO D o" CD | i a —-s D o a 0 o Q 8 CV CD ^ CO o o o O o ^ !§ c 3 cr CD CD O 3 a CD <D 3 3 o 0) it* CD Q. Committee's Page. This page is not to be hand-written except for the signatures 1 To my parents Qin and Zhao. 11 ACKNOWLEDGMENTS It was quite a journey pursuing the Ph.D., and I would not have been able to complete this journey without the support of many people over the past several years. I must first express my gratitude towards my advisors, Professor Susan Hagness and Professor Daniel van der Weide. They not only helped me keep growing technically and professionally, but their leadership, vision, and strategic thinking also set an example I hope to match some day. I wish to thank those faculties and technical staff at UW for their consulting and assistance in technical areas: Alan Bettermann, Robert Blick, Timothy Hall, Steve Limbach, Robert Marsland, Charles Paulson, Barry Van Veen, and Tomy Varghese. I would like to thank the graduate students I have worked with on different projects in UWCEM and vdW groups: Min-Ki Choi, Anuj Dron, Hong-Joon Kim, Al Mashal, and Jake Shea. I look forward to future collaboration with any of them. I also thank some of my fellow Ph.D. students: Suzette Aguilar, Matthew Burfemdt, Zhen Ji, Steve Kennedy, Mariya Lazebnik, Keely Willis, and Earl Zastrow. They each helped make my time in the Ph.D. program more fun and interesting. Last but not least, millions of thanks to my parents Qin and Zhao who have been supporting every decision I made over the years. iii TABLE OF CONTENTS Page LIST OF TABLES v LIST OF FIGURES vi ABSTRACT 1 2 1 1.1 1.2 1.3 1.4 1.5 1 2 3 4 5 Background of Current Breast Cancer Detection Techniques Overview of Active Microwave Imaging Techniques for Breast Cancer Detection . Motivation to Develop a Low-Cost Microwave Refiectometry Motivation to Investigate an Acoustic and Microwave Hybrid Modality . Contributions Ultrawideband Compact Mixer-based Reflctometer (CMR) Background Knowledge Theory and Design Calibration Fabrication and Measurement Results . Motivation and Background Free-Space CMR-Antenna Measurement System Free-Space Calibration Method Results Concept of Acoustic and Microwave Hybrid Modality 4.1 4.2 7 . Calibrated Free-Space Measurements with A CMR-Antenna Subsystem 3.1 3.2 3.3 3.4 4 xi Introduction 2.1 2.2 2.3 2.4 2.5 3 .:.... Theory . . Methods of Inducing Tissue Deformation 7 8 10 11 12 15 15 16 16 20 26 26 27 IV Page 5 Global-Acoustic and Microwave Hybrid Modality 5.1 5.2 5.3 5.4 6 7 8 9 2-D Numerical Simulation Discussion of Limitations 3-D Numerical Simulation with a Bulk Translation Removal Method Discussion of Limitations 29 . 29 37 37 46 Focused-Acoustic and Microwave Hybrid Modality 6.1 Local Excitation with a Focused Acoustic Transducer 6.2 Overview of the 3-D Numerical Techniques 49 . 49 49 3-D Simulation of Focused-Acoustic Induced Tissue Deformation 52 7.1 Acoustic Beamformer Design and Modeling 52 7.2 Induced Deformation Results 54 3-D Simulation of Microwave Interaction with a Deforming Tissue Target 8.1 Derivation of 3-D SBCs 8.2 3-D FDTD-SBC Implementation 8.3 Validation 3-D Computational Investigation of Focused-Acoustic and Microwave Hybrid Modality ...... 9.1 Results - Microwave Scattering from a Canonically Deforming Target 9.2 Results - Microwave Scattering from a Focused-Acoustic Induced Deforming Target 9.3 Results - Impact of Dielectric and Elastic Properties Variations on the Microwave Doppler Scattering -....' 60 60 66 68 74 74 75 78 10 Conclusion and Future Work 83 LIST OF REFERENCES . . . . . . . 86 APPENDICES Appendix A: Appendix B: Appendix C: Theoretical Analysis of the Doppler Reconstruction Method List of Authored Publications Developed Codes for FDTD-SBC Implementation 93 96 98 V LIST OF TABLES Table 5.1 Page Introduced errors in calculating the bulk translation / 47 VI LIST OF FIGURES Figure Page 2.1 Signal flow analysis of a reflectometer 7 2.2 Mixer-based 6dB bridge-T attenuator circuit 8 2.3 Mixer-based reflectometer system with two external signal sources 9 2.4 A prototype of the cascaded 6dB bridge-T attenuators 2.5 (a) Magnitude and (b) Phase of the reflection coefficients of the DUTs measured by the reflectometer prototype as well as by a commercial PNA 13 2.6 Simulated and measured directivity and coupling factor for the reflectometer prototype 14 3.1 CMR-antenna system: The compact CMR and a UWB horn antenna are connected through a semi-rigid cable 16 Prototype: The compact CMR and a UWB horn antenna are connected through a semi-rigid cable 17 3.2 . . 10 3.3 The controlled experiment setup used in the free-space calibration and measurement procedure . 19 3.4 Raw time-domain responses from the standards and MUT 21 3.5 Time-gated and calibrated time-domain responses from the standards and MUT . . . . 22 3.6 Measured and calibrated magnitude of complex reflection coefficients of a silicone rubber sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMR-antenna system are compared with those obtained using commercial PNA-antenna setup as well as with the theoretical values 22 Vll Figure 3.7 Page Measured and calibrated phase of complex reflection coefficients of a silicone rubber sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMRantenna system are compared with those obtained using commercial PNA-antenna setup as well as with the theoretical values. 23 Measured reflection coefficients of (a) a Silicone rubber sheet (er = 3, thickness = 1.2cm) and (b) a Skycast material (er = 12, thickness = 2.6cm) using the free-space CMR-antenna prototype. Results (circles) compare favorably with the theory values (dotted lines) 24 Measured reflection coefficients of (c) an FR4 board (er = 4.34, thickness = 0.2cm) and (d) a Plexiglas material (er = 2.7, thickness = 2.4cm) using the free-space CMRantenna prototype. Results (circles) compare favorably with the theoretical values (dotted lines) 25 Potential ways of inducing tissue deformation in breast: (a) global excitation with a compression plate, and (b) local excitation with a focused acoustic transducer 28 Comparison of the amplitudes of the normalized Doppler component computed via FDTD-SBCs (dashed line) and the normalized difference between stationary scattered field solutions computed via FDTD-REP (solid line). The 2D scattering object is a deforming 4-mm-diameter circular cylinder undergoing a peak radial displacement of 0.45 mm 33 Phantom configuration for the hybrid simulations, with axis of cylindrical inclusion into the page . 33 Deformed inclusion boundaries, under a 10% peak applied strain (Bulk translation is removed). Peak deformation of boundary point (0, 2) is 0.24 mm for the normal fibroglandular inclusion and 0.07 mm for the malignant tumor inclusion 35 Computed Doppler components in the backward (0 = 180°) and forward (</> = 0°) scattered fields under different peak applied strains, for malignant and normal fibroglandular scatterers 36 3-D numerical phantom testbed for the (a) mechanical simulation and (b) electromagnetic simulation of the global-acoustic and microwave hybrid modality 39 5.6 An illustration of the 'null' detection procedure 41 5.7 An illustration of the bulk translation calculation procedure 41 3.8 3.9 4.1 5.1 5.2 5.3 5.4 5.5 Vlll Appendix Figure 5.8 5.9 Calibrated elastic boundary deformation of the malignant tumor and the normal fibroglandular tissue clutter. Bulk translation for both cases are removed with a synthetic tracking method through the use of an antenna array. The normal fibroglandular tissue clutter experiences about three times the boundary deformation of the malignant tumor. 7.1 7.2 43 The Doppler scattering response from the deforming malignant tumor and deforming normal fibroglandular tissues. Despite the dielectric properties variation in the normal fibroglandular tissue, the Doppler scattered field from the deformed malignant tumor is observed to be 5-10 dB weaker than that from the deformed normal fibroglandular tissue clutter . 44 5.10 The Doppler scattering response from a deforming malignant tumor and a deforming normal fibroglandular tissue clutter, when different levels of errors in approximating the bulk translation are introduced in the testbed 6.1 Page 45 Computational steps of simulating microwave scattering response from a deforming target under focused-acoustic excitation 51 Acoustic intensity (W/cm2) in the (a) lateral and (b) axial focal plane generated by a focused square transducer in a normal fibroglandular tissue background: The achieved focal region is roughly 1 cm x 1 cm x 3.5 cm 55 3-D computational testbed for simulating the microwave Doppler scattering response from a tissue target excited by a focused acoustic transducer. The target is either a malignant tumor or a normal fibroglandular tissue clutter, surrounded by normal breast tissues 56 7.3 Focused-acoustic induced time-varying radial deformation of a tumor boundary, ipiaterai is the angle with respect to the y-axis in the y-z plane and ipaxiai is the angle with respect to the x-axis in the x-z plane. The group of solid-line curves represents deformation for different snapshots on the acoustic time scale. Its largest boundary deformation is 0.07 mm 57 7.4 Focused-acoustic induced time-varying radial deformations of a normal fibroglandular tissue clutter boundary, ipiaterai is the angle with respect to the y-axis in the y-z plane and ipaxiai is the angle with respect to the x-axis in the x-z plane. The group of solidline curves represents deformation for different snapshots on the acoustic time scale. Its largest boundary deformation is 0.23 mm 58 ix Appendix Figure 8.1 8.2 8.3 8.4 8.5 8.6 8.7 9.1 -.. Page Sheet impedances Zs and admittances Ys for a 3-D dielectric object are derived by enforcing equivalence in the scattered field solutions for a) a perturbed (r = a + 5) dielectric sphere, and b) an unperturbed (r = a) dielectric sphere with SBCs applied at the boundaries. 61 3-D numerical testbed for validating the FDTD-SBC implementation: simulate the microwave scattering response from a deforming malignant tumor surrounded by a normal fibroglandular tissue background 69 The computed total scattered fields from a slightly bloated spherical tumor with a uniform radial perturbation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, ip = —90° is the monostatic angle 69 The computed total scattered fields from a perturbed spherical tumor with a nonuniform radial perturbation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, ip — —90° is the monostatic angle. 70 The computed total scattered fields from a perturbed ellipsoidal tumor with a nonuniform radial perturbation across its boundary, -p is the angle with respect to the y-axis in the y-z plane. tp_ = —90° is the monostatic angle 70 The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor with a uniform peak radial deformation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, tp = —90° is the monostatic angle 72 The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor with a non-uniform peak radial deformation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, ip = —90° is the monostatic angle 73 Simulated magnitude of Doppler scattering responses of a canonically deforming target, tpiaterai is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90° is the monostatic angle. There is a 6-8 dB contrast in the Doppler signal strength between the uniformly-deforming malignant tumor (5=0.1 mm) and normal fibroglandular tissue clutter (5=0.3 mm) cases. This contrast is larger and inverted, in comparison with the contrast in fundamental scattering from the two types of targets as in pure microwave detection techniques 76 gure Pag 2 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced deforming target, ipiaterai is the angle with respect to the y-axis in the y-z plane and ^lateral = —90° is the monostatic angle. The Is* Doppler scattering contrast between the deforming malignant tumor (5max=Q.01 mm) and normal fibroglandular tissue clutter (5max=Q.2?> mm) cases is larger and inverted, in comparison with the contrast in the fundamental scattering from the two types of targets as in pure microwave detection techniques. 77 3 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced deforming target which varies in its dielectric properties. The assigned dielectric properties are: the 50th percentile data of normal fibroglandular tissue (e r -44, CT=4.16), the 75th percentile data of normal fibroglandular tissue (er=48, (7=4.71), and the 50th percentile data of malignant tumor (er=50, <r=4.91). ^lateral is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90° is the monostatic angle . . . 81 4 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced deforming target which varies in its elastic properties. The assigned elastic properties are: E=12 kPa, 9 kPa, and 7.5 kPa that represents different types of malignant tumors, and E=1.82 kPa that represents a normal fibroglandular tissue clutter, v = 0.495 is assigned to all these target cases, ipiaterai is the angle with respect to the y-axis in the y-z plane and ^; atera ; =—90° is the monostatic angle. . . 82 XI ABSTRACT This dissertation focuses on technology development that contributes to cost reduction and sensitivity improvement in microwave breast cancer detection systems. A microwave refiectometer measures scattered signals from breast tissue, which allows subsequent signal processing to obtain breast images for screening and diagnosis purposes. In current systems, a commercial vector network analyzer (VNA) functions as a refiectometer, which connects to either a single antenna or an antenna array via switches. An ideal configuration would be a fully-populated reflectometer-antenna array with each antenna connected to its own refiectometer. The VNA is too bulky and expensive for this approach. To address this problem, a low-cost mixerbased reflectometer-antenna subsystem is developed. A free-space calibration method is developed to eliminate the need for connectorized mechanical standards, compensate for antenna reverberation, and allow flexibility in the choice of reference plane position in measurements. Measurement results using the prototype compare favorably with those using a VNA. This free-space calibrated reflectometer-antenna subsystem offers potential for cost reduction in a full-array based microwave breast imaging system. Heterogeneity in the dielectric properties of normal breast tissue poses a challenge for current microwave techniques to distinguish malignant from normal fibroglandular tissue. Hence, an acoustic and microwave hybrid modality that exploits both dielectric and elastic properties contrasts in breast tissue is proposed. Two methods of inducing tissue deformation in breast — global excitation with a compression plate and local excitation with a focused acoustic transducer — are Xll investigated. 3-D numerical techniques are developed to efficiently model this multi-physics problem. An acoustic beamformer is designed and the finite-element-method based acoustic-structural simulation is conducted to compute the induced tissue deformation. The finite-difference timedomain based numerical method with the use of sheet boundary conditions is adopted and extended to 3-D to simulate the microwave scattering response of the acoustically excited tissue. The numerical study demonstrates a larger and inverted microwave Doppler scattering contrast between the malignant and normal fibroglandular tissues in comparison with their fundamental scattering contrast of a pure microwave detection scheme. This reveals the potential of the hybrid modality for improved sensitivity in detecting a malignant tumor in a heterogeneously dense breast environment. 1 Chapter 1 Introduction 1.1 Background of Current Breast Cancer Detection Techniques Breast cancer is the second most common cause of death from cancer and the most common non-skin-related malignancy among women in the United States, according to the US Center of Disease Control and Prevention. In 2005, over 180,000 new female and over 1,700 new male cases of invasive breast cancer were diagnosed; in the same year, more than 40,000 women and more than 370 men died from the disease. Early detection of malignancy, based on regular screening and diagnosis programs, and thus timely medical treatment is very important in saving lives of the breast cancer patients. X-ray breast mammography is currently the primary screening tool for detecting non-palpable early-stage breast cancer, and is recommended for women aged 35 and older in the United States. Technical advances such as digital mammography, as well as advances in radiological expertise, have greatly improved the overall quality and interpretation of mammographic images. However, the sensitivity-related limitations of mammography result in a relatively high false negative rate (4-34%) [1], particularly in patients with radiographically dense breast tissue such as young women [2]. Furthermore, the relatively high false positive rate (70%) [3] of mammography leads to many unnecessary biopsies. Other drawbacks of X-ray mammography include uncomfortable or painful breast compression and exposure to low levels of ionizing radiation. Ultrasound scanning, also called sonography, is one of the main alternative imaging techniques used in clinics [4, 5]. The ultrasound breast exam does not use ionizing radiation, is non-invasive, and captures images in real-time. The main use of the ultrasound procedure is determining if an 2 abnormality is solid (which may be a non-cancerous lump of tissue or a cancerous tumor) or fluidfilled (such as a benign cyst) or both cystic and solid. However, ultrasound scanning alone cannot determine the nature of a solid abnormality. Contrast-enhanced MRI is another main alternative imaging technique for breast cancer detection and is playing an important role in the diagnostic evaluation of mammographically detected lesions [4, 5]. Contrast-enhanced MRI is minimally-invasive, has good resolution, and provides valuable information for diagnostic purposes. Yet, this modality is quite costly for routine screening purposes, takes longer to perform than other modalities, and is not recommended for patients with metallic medical implants. The tremendous toll that breast cancer takes, combined with the limitations of the current methods, has motivated intense research efforts to advance alternative techniques as complimentary diagnostic tools. 1.2 Overview of Active Microwave Imaging Techniques for Breast Cancer Detection Active microwave techniques that exploit the dielectric-properties contrast in breast tissue have gained increased research interests in the past decade for breast cancer detection. The fundamental physics that serves as the basis for the active microwave techniques is the dielectric-properties contrast in breast tissue in the microwave frequency range suggested by several dielectric spectroscopy studies over the past several decades [6, 7, 8, 9, 10]. In microwave breast imaging techniques, tissue-dependent microwave scattered signals, received by antennas and measured by microwave refiectometers, are processed to obtain breast images for screening and diagnosis purposes. Microwave tomography [11, 12, 13, 14, 15, 16] and microwave radar [17, 18, 19, 20, 21, 22] are the two main classes of current active microwave imaging techniques. Microwave tomography recovers the dielectric-properties profile of an object from measurement of the transmission of microwave energy through the object. This requires the solution of an ill-conditioned nonlinear inverse-scattering problem via elaborate image reconstruction algorithms. Some promising initial results for breast imaging have been obtained under specialized 3 conditions in 2-D [13], and recent advances in 3-D microwave inverse scattering have been reported in literature [14, 15, 16]. Microwave radar imaging technique uses scattered signals to infer the locations of significant scatterers in breast. Significant scattering arises from significant contrasts in the dielectric properties, for example, between normal adipose dominated tissues and malignant tumors. This modality departures from the complicated image reconstruction techniques required by conventional microwave tomography, as it only seeks to identify the presence and location of strong scatterers in the breast, rather than to completely reconstruct the dielectric-properties profile. Several promising experimental studies have been reported [21, 22]. Both the active microwave imaging techniques exploit the dielectric properties contrast in breast tissue, and the current prototypes in both techniques employ a commercial VNA to function as the refiectometer for measuring microwave scattered signals [13, 21, 22]. 1.3 Motivation to Develop a Low-Cost Microwave Reflectometry A microwave refiectometer, a key hardware in microwave breast cancer detection systems, can measure the tissue-dependent microwave scattered signals, which allows subsequent signal processing to obtain breast images for screening, diagnosis and monitoring purposes. In current systems, one single commercial VNA is employed to function as a high performance refiectometer, which connects either to a mechanically-scanned single antenna or to an antenna array via switches. However, both approaches have their limitations. The effect of cable movements during the antenna scanning process is found to be a significant source of measurement non-repeatability. For the latter approach, additional calibration errors would be introduced, and the switching implementation would limit the employment of multi-channel transmit-receive strategies that offer potential improvement in resolution and robustness to artifacts in the imaging systems. An ideal configuration instead would be a fully-populated refiectometer-antenna array with each antenna connected to its own refiectometer. VNA is too bulky and expensive for this approach, which 4 motivates the development of a low-cost and compact reflectometer that can be paired with an ultrawideband (UWB) antenna as a rigid subsystem for receiving and measuring microwave scattered signals. Different methods were proposed to implement a microwave reflectometer. A low-cost reflectometer was constructed using planer directional couplers which are composed of cascaded Wilkinson power dividers [23]. However, this is a narrowband design and therefore is not suitable for an ultrawideband application. A six-port architecture based on power detection was proposed for microwave network analysis, but it suffers from the problem of low dynamic range which is a side effect of employing power detection instead of coherent detection, limited bandwidth and difficulty of calibration [24, 25, 26, 27]. Wheatstone bridge was also used for reflection coefficient measurement and has its application in mm-wave instrumentation area [28, 29]. 1.4 Motivation to Investigate an Acoustic and Microwave Hybrid Modality As mentioned in Section 1.2, the fundamental physics that serves as the basis for active mi- crowave techniques is the dielectric-properties contrast in breast tissue in the microwave frequency range. Among those reported dielectric spectroscopy studies [6, 7, 8, 9, 10], early studies [6, 7, 8] characterized the dielectric properties for a general group of normal tissue, whereas more recent studies [9, 10] broke down the normal tissue characterization into groups with different adipose content. It is known from the studies that the dielectric properties contrast between malignant and normal adipose dominated tissues is as large as on the order of 10:1; However, the contrast between malignant and normal fibroglandular tissues is small, about 10%. A human breast is a dielectrically-heterogeneous tissue environment, and very often tumors grow in or near normal fibroglandular tissues. Therefore, researchers are seeking novel ways to enhance the contrast between malignant and normal fibroglandular tissues in order to achieve better detection sensitivity. In the past decades, researchers have also conducted experimental studies to measure the elastic modulus of human breast tissue. Reported measurements from different research groups were often conducted following different procedures which do not necessarily produce exactly the same results on the same type of tissue. In addition, different tissue handling techniques have been 5 used in generating the data, which is also a source of uncertainty in evaluating the results found in literature. But overall, despite the above variabilities, previous studies [30, 31, 32, 33, 34] indicate a significant elasticity contrast between malignant tumor and normal breast tissues (both adipose-dominated and fibroglandular types). As a brief summary, the Young's modulus of different carcinomas are found to be much larger than that of either normal fibroglandular tissue or normal adipose tissue [30, 31]. Shear modulus and shear viscosity in breast cancer are found to be higher than those in surrounding normal breast tissues [32, 33]. Invasive and infiltrating ductal carcinoma is found to be significantly stirrer than normal glandular tissues [34]. Such elastic contrast has motivated research effeorts in various elastography applications for breast imaging, including breast sonoelastography [35, 36], breast magnetic resonance elastography (MRE) [37], and breast vibro-acoustography [38]. The author is thus motivated to investigate an acoustic and microwave hybrid modality that exploits both dielectric and elastic properties contrasts in breast tissue. Microwave-induced thermoacoustic tomography was proposed for breast cancer detection as in [39,40,41,42]. However, the physics basis of the microwave-induced thermoacoustic tomography is very different from the proposed acoustic and microwave hybrid approach presented in this dissertation. In microwave-induced thermoacoustic tomography, microwave energy is selectively absorbed by the high conductivity tissues in breast and induces acoustic waves by means of thermoelastic expansion. Ultrasound transducers are used to detect the thermoacoustic signals which are then processed for image reconstruction. 1.5 Contributions In this dissertation, a low-cost mixer-based reflectometer-antenna subsystem is developed. The novel mixer-based reflectometer design allows significant cost and size reduction. A free-space calibration method is developed to improve the measurement repeatability by eliminating the need for contacting mechanical standards, compensating for antenna reverberation, and allowing flexibility in the choice of reference plane positions in measurements. Measurement results using the prototype compare favorably with those using a VNA. This free-space calibrated reflectometer-antenna 6 subsystem offers potential for cost reduction in a full-array based microwave breast imaging system. In this dissertation, an acoustic and microwave hybrid modality that exploits both dielectric and elastic properties contrasts in breast tissue is also proposed. Two methods of inducing tissue deformation in breast — global excitation with a compression plate and local excitation with a focused acoustic transducer — are investigated. 3-D numerical techniques are developed to efficiently model this multi-physics problem. An acoustic beamformer is designed and finite-elementmethod (FEM) based acoustic-structural simulations are conducted to compute the induced tissue deformation. The finite-difference time-domain (FDTD) based numerical method with the use of sheet boundary conditions (SBCs) is adopted and extended to 3-D to simulate the microwave scattering response of the acoustically excited tissue. From our numerical studies, we found that the local excitation approach has its limitations in order to be used in this hybrid modality for breast cancer detection such as the presence of the bulk translation and the difficulty of removing the bulk translation with a practical antenna array. The focused-acoustic and microwave hybrid modality, however, overcomes this problem. The numerical simulation demonstrates a larger and inverted microwave Doppler scattering contrast between the malignant and normal fibroglandular tissues in comparison with their fundamental scattering contrast of a pure microwave detection scheme. This reveals the potential of the hybrid modality for improved sensitivity in detecting a malignant tumor in a heterogeneously dense breast environment. 7 Chapter 2 Ultrawideband Compact Mixer-based Reflctometer (CMR) This chapter presents the design, fabrication, calibration and measurement results of an ultrawideband compact mixer-based reflectometer (CMR). 2.1 Background Knowledge \ k L D-directivity C-coupling 3 2 UUT 2 Figure 2.1 Signal flow analysis of a reflectometer A reflectometer is a circuit that uses a wave separator to isolate and sample the incident and reflected waves from a mismatched load. It measures the reflection coefficient of a one-port network and, in a more general configuration, the S parameter of a two-port network. Figure 2.1 shows the signal flow analysis of a reflectometer [43]. In operation, the separator provides a sample Vi of the incident wave and a sample Vr of the reflected wave. A realistic reflectometer has finite directivity D, which means both the incident wave and the reflected wave contribute to both Vi and Vr, and therefore an error in the reflection coefficient measurement is introduced. Assume a unit incident wave, the expressions of Vi and Vr can be obtained from the signal flow analysis as Vi = C C, 9 TeP , D C_ + CTe?* D Vr (2.1) where T is the complex reflection coefficient of the device under test (DUT). 9 and 4> are the phase delay difference through the reflectometer circuit, C and D are the coupling factor and the directivity of the reflectometer respectively. The calibrated ratio Vr/Vi can represent the reflection coefficient of the DUT as I v. I I r l +4 (2.2). V The measurement uncertainty due to finite directivity is approximately ± (1+ | F |) /D, which implies that higher directivity results in better measurement accuracy. 2.2 Theory and Design R1 Zin ZS *«^/\^u~~» R25 Zout: -WW < ^ ^ > Figure 2.2 Mixer-based 6dB bridge-T attenuator circuit A resistive attenuator circuit is employed as a single wave separator in this reflectometer design to achieve compactness. Based on the simulation using Agilent Advanced Design System (ADS) software, a 6 dB bridge-T attenuator is chosen because it maintains good port characteristic impedance over the 1-12 GHz frequency range used in this ultrawideband microwave imaging application. Furthermore, in a novel arrangement, the normal 50 Q, resistors to ground are replaced by frequency mixers with a similar input impedance instead as shown in Figure 2.2. This reflectometer enables direct detection and measurement of the complex reflection coefficients while absorbing 9 the parasitic resistance of the active elements into part of the attenuator circuit to dramatically reduce the effects of parasitic. Because the broadband mixers chosen for construction do not have a characteristic impedance of Z0 at the RF input port, a 3 dB attenuator is placed between the resistive part and the RF input port of the mixer to increase isolation for better broadband performance. VvV •AM- OUT (V,RF 3dBAttn 3dB Attn Computer { N ' IF1 "lo \'S} IF2 Reference 1 Loefeirv. Amplifier Figure 2.3 Mixer-based reflectometer system with two external signal sources Two cascaded bridge-T attenuator circuits function to separate and detect samples of the incident and reflected waves. In this coherent detection system as shown in Figure 2.3, two synchronized external sources provide RF and LO signals and a lock-in-amplifier measures the complex IF signals detected by the mixers. The expressions for the two IF signals can be obtained quantitatively through microwave theory analysis as [-V + -Vrei29}Aei2e\ IF1 = IF2 = [^Vejd + -VTe^9]Bej29" 4 4 (2.3) where V represents the input RF signal, T is the complex reflection coefficient of the DUT, 9 is the electrical length of the transmission line connected between two attenuator circuits, Ae^26' and Bej2e" represent the contribution of LO signal to IF output signals. The signals IF1 and IF2 both depend on the incident wave as well as the reflected wave; however, IF1 depends relatively more on the incident wave whereas IF2 depends more on the reflected wave. Thus, the difference between the two IF signals after calibration can represent the complex reflection coefficient T: IF1-K-IF2 oc T (2.4) 10 where K is the complex calibration factor. This complex reflection coefficient is referred as the raw reflection coefficient Traw in the following discussions. 2.3 Calibration A Figure 2.4 A prototype of the cascaded 6dB bridge-T attenuators Three types of measurement errors are introduced in the microwave refiectometer system: systematic errors, random errors, and drift errors. Systematic errors are caused by imperfection in the refiectometer setup, for example, the imperfect broadband performance of the bridge-T attenuators, the mismatches at the RF input port, DUT port and IF output ports, the reflection between the resistive circuit and the input port of frequency mixer, and the imperfection due to multiple use of adapters and connectors, etc. If these systematic errors do not vary over time, assuming that a one-time measurement is short enough in time, then the systematic errors can be removed mathematically through a calibration procedure. Random errors vary as a function of time and cannot be removed mathematically since they are unpredictable. In the refiectometer system, the main sources of random errors are the instrument noise (phase noise of the source signal and sample noise in the lock-in-amplifier) and the connector n non-repeatability. Random errors yet can be reduced by increasing the input source power, narrowing the system IF bandwidth and averaging over multiple sweeps. Trial experiment with the reflectometer prototype has been conducted to find the optimal driven LO and RF power as well as the appropriate IF bandwidth and averaging without sacrificing too much measurement time. Drift errors occur when the system performance changes after a calibration has been performed, which is mainly caused by temperature variation and mechanical instability. The thermal characteristics of multiple cables used in the reflectometer system make the system sensitive to drift errors, especially when the system becomes heated after being in operation for hours. Efforts are made to ensure the mechanical stability of the whole experimental setup to further reduce the drift errors. A one-port error network models the systematic errors which are categorized into three main terms: ED, ERT, and Es. ED is the directivity error due to signal leakage, ERT is the frequency response error caused by reflection and transmission tracking within the system, and Es is the error due to source or impedance mismatches. These three error terms can be derived from the general equation below by establishing three equations with three known DUTs: TA = ED + ERT • F X ™ III (2.5) 1 raw To establish the equations, three known calibration standards need to be measured. In the reflectometer system, an open, a short and a load (50 fi), which is known as the SOL calibration standards, are used. Solving these equations then yields the systematic error terms and makes it possible to derive the actual reflection coefficient r ^ of a DUT from measured raw reflection coefficient Traw. 2.4 Fabrication and Measurement A mixer-based reflectometer measurement system is shown in Figure 2.3. The cascaded 6 dB bridge-T attenuators as shown in Figure 2.4 are fabricated on a Rogers RO3010 printed circuit board (PCB) (er = 10.2) with surface mount chip resistors. This microstrip line design has an approximately 2 cm x 2 cm planar size and thus can accommodate the size of SMA connectors. A Marki Ml-0212 frequency mixer which covers the RF frequency range 2-12 GHz is selected for 12 frequency conversion. An SRS844 RF lock-in-amplifier receives and detects the complex value of IF1 and IF2 signals. For this low power detection application, a +9 dbm RF input signal is chosen to feed the reflectometer which, after attenuation through the cascaded 6 dB bridge-T circuits, produces an approximately -3 dbm power signal ready to be transmitted to the DUT. LO driven power for the mixer is chosen to be +16 dbm for the optimal linearity performance. Mechanical open, short and load standards from Agilent 85033D calibration kit are used to perform the oneport SOL calibration. Two different DUTs are used for measurement: one is an open-ended 3 dB attenuator, and the other is a DUT composed of several connectors and fixed-value termination. Since the exact reflection coefficient of the second DUT is not known intuitively, it is referred as an 'unknown DUT' in the following discussions. 2.5 Results The measured complex reflection coefficient results for the two DUTs are summarized in Figure 2.5. The measurement results using the reflectometer prototype show a fairly good agreement with the results measured by a commercial Agilent E8364A Performance Network Analyzer (PNA). In addition, directivity and coupling factor, as the important performance indices of a reflectometer, has also been simulated in ADS and measured in the lab. As shown in Figure 2.6, the coupling factor of the reflectometer is 7 dB ± 1.5 dB, and its directivity is around 6-12 dB in the wide frequency range of 1-12 GHz. In the next chapter, this developed compact mixer-based reflectometer technique is integrated into an ultrawideband CMR-antenna subsystem for free space microwave measurements. 13 -5 vw3?£f&*££?&>' §-10 -15 -20 o 3dB atteunator DUT (Our Prototype) « Unknown DUT (Our Prototype) 3dB attenuator DUT (PNA) Unknown DUT (PNA) 4 6 8 Frequency (GHz) 10 12 10 12 (a) 4 6 8 Frequency (GHz) (b) Figure 2.5 (a) Magnitude and (b) Phase of the reflection coefficients of the DUTs measured by the reflectometer prototype as well as by a commercial PNA 14 16 \A Directivity Measurement —«— Directivity Simulation Coupling factor Measurement Coupling factor Simulation 121S10; 6 8 frequency (GHz) Figure 2.6 Simulated and measured directivity and coupling factor for the refiectometer prototype 15 Chapter 3 Calibrated Free-Space Measurements with A CMR-Antenna Subsystem This chapter presents a free-space calibration method and the free-space calibrated microwave measurements using an ultrawideband CMR-antenna subsystem. 3.1 Motivation and Background Current microwave imaging systems under development typically employ antennas in combination with a single reflectometer, mostly a commercial VNA, to receive and measure scattered signals [13, 21, 22]. Traditional calibration procedures used in these systems require repeated mechanical contact and do not compensate for antenna dispersion, mechanical artifacts or reverberation, which motivates work to investigate a free-space calibration method that is suitable for the CMR-antenna subsystem in order to avoid the above problems. A variety of free space measurement systems and calibration techniques have been reported for use in the microwave regime [44, 45, 46, 47, 48]. The measurement systems in [44, 45] are based on a two-port thru-reflect-line (TRL) calibration. A modified reduced-cost measurement system which follows this approach is reported in [46]. Reference [47] constructs a two-parameter model to approximate antenna-media coupling and uses shifted metal plates as calibration standards. Reference [48] discusses free space calibration procedures using two-term and three-term error models, as well as free space calibration standards and errors associated with the calibration methods. 16 3.2 Free-Space CMR-Antenna Measurement System ref Lock-in PC t IF1 ^ S f-C^Zyd JVW Figure 3.1 CMR-antenna system: The compact CMR and a UWB horn antenna are connected through a semi-rigid cable. A calibrated free-space measurement procedure is developed for a monostatic system comprised of a microwave CMR and a UWB horn antenna. The calibration procedure involves a precalibration using the empty room standard case, a time-gating process, and a three-term one-port error model. Multiple versions of the low-cost CMR circuit block, as shown in Figure 3.1, can be arranged in parallel for array-based systems. This contactless measurement system eliminates the need for connectorized mechanical standards at the reference planes, thus providing better measurement reliability and repeatability. The calibration method compensates for the UWB antenna dispersion, the reverberation from the antenna, and reflections between the antenna and the material under test (MUT). It also allows for flexibility in the choice of reference plane positions. Figure 3.2 shows a combined CMR and UWB horn antenna prototype. 3.3 Free-Space Calibration Method The measurement and calibration procedure used to compensate for the error sources and measurement artifacts inherent in the whole system is described below. First, extract the broadband pre-calibrated measured reflection coefficient Traw from the raw outputs IF1 and IF2 of the CMR. X is a pre-calibration coefficient determined by forcing the Fraw of the empty room case to be zero. rP„,„ = X-IF1- IF2 (3.1) 17 ^ I I Figure 3.2 Prototype: The compact CMR and a UWB horn antenna are connected through a semi-rigid cable. 18 Second, apply a time-gating step to obtain a smoothed reflection coefficient FQ- Only the timedomain reflection contributions from the standards or MUT is of interest, thus a window function H(t) is selected which excludes any other reverberation or background noise and clutter that fall outside the window. This window is applied to the inverse discrete Fourier transformed timedomain data ^ r _ 1 {r r . o l u }. As a result, a time-gated response TG is obtained by Fourier transforming the time-domain gated information back into the frequency domain as YG = T{T-l{Yraw}-H{t)} (3.2) Finally, apply a one-port error model to extract the free space system error terms ei, e2, e3, and to obtain the calibrated reflection coefficient Tc- rc = - Tc-ei (3.3) e2 + e 3 • ( r G - ei) In the free-space measurement system, e\ represents the return loss present in the CMR circuit and the mismatch between the CMR and antenna; e2 accounts for any UWB antenna dispersion, multipath effects, semi-rigid cable loss and free space path loss; e 3 takes into account the multiple reflections returning from the MUT. To find the error terms, three known calibration standards are employed: 'Empty room' (load), a flat copper sheet placed at the reference plane (reference short), and the same copper sheet shifted a distance Ad away from the reference plane (offset short), as illustrated in Figure 3.3. The reference plane is placed in the far field region for this application. Finally, the three error terms are obtained and used to determine the unknown material reflection coefficient calibration, TER, TRS> and Tos are the time gated reflection coefficients of the three standards. 1 0 0 ei 0 — {FRS — FER) -1 e2 FRS— e3 TQS 0 _e-3P2Ad{Tos __ TER) _e-jf32Ad FER _ FER FER (3.4) 19 y y y y y y y v j v yV y v y v v a) CMR K] S < < < i i < < ^AiVsAiiiiai^ b) £sSn*!\£iA 1 Copper sheet<] Cft/IR yMi^iUi^Sw^^^^ 7^?V?;7"^ r-.^v'v^TVVvvf.v-^ M c) CIViR /vVSfTv^ C) CIMR Reference plane Figure 3.3 The controlled experiment setup used in the free-space calibration and measurement procedure 20 3.4 Results Compact (2 cm x 2 cm) microstrip-line based bridge-T attenuator circuits as described in Chapter 2 were fabricated, and a UWB horn antenna was connected to the output port of circuits through a 30.5 cm semi-rigid SMA cable. The antenna used was the Tera-X Model 530, a 4.2 cm x 2.6 cm modified horn antenna that is based in part on work in [49]. The RF input source was swept from 1 to 12 GHz with a step size of 50 MHz, yielding 221 data points. The RF power was optimized at +8 dbm for better wideband performance. As shown in Figure 3.3, the free-space calibration and measurement area was shielded with absorbing foam (ECCOSORB-AN) to provide an anechoic environment (61 cm x 76.2 cm x 76.2 cm). Optical rails were also used to achieve precisely controlled and repeatable standard placement. The desirable reference plane position is in the far-field region for measurement. However, the reference plane position needs to be closer enough to the transmit antenna in order to achieve a reasonable signal-to-noise (SNR) ratio. With this trade-off considered, we chose the reference plane to be 10 cm away from the antenna aperture. A 30.5 cm x 30.5 cm copper sheet was used as the reference short and the offset short standards; whereas a 30.5 cm x 30.5 cm flat silicone rubber sheet (er = 3, thickness = 1.2 cm, Durometer = 70 A) served as the MUT. The calibrated time-domain responses in Figure 3.4 and Figure 3.5 show a vivid picture of the effectiveness of the free-space calibration. The inverse discrete Fourier transformed time-domain return pulses from the standards as well as the MUT, both prior to calibration as shown in Figure 3.4, i.e., .F -1 {rV ot0 }, and after time-gating and free-space calibration as shown in Figure 3.5, i.e., Jr~1{Tc} are shown. A bandpass Gaussian window is used in the transformation. Without calibration, the original reference plane for the time-domain observation is at the output port of the CMR. The pulses complete a round trip, traveling through the semi-rigid cable and over the free-space path, then being reflected and returning. As in Figure 3.4, they arrive approximately 4.3 ns after passing the CMR reference plane. In Figure 3.5, the calibrated return pulses are observed 21 immediately after a time delay defined for the modulated Gaussian source pulse because the reference plane has been relocated by the free-space calibration procedure. It is apparent that the pulse broadening due to dispersion in the circuit and antenna is effectively compensated. 1] 5 Time (ns) Figure 3.4 Raw time-domain responses from the standards and MUT Figure 3.6 and Figure 3.7 show the calibrated complex reflection coefficients of a silicone rubber sheet measured using the CMR-antenna system. This result compares favorably with those obtained using an Agilent E8364A PNA-antenna setup in the frequency range of 3 to 12 GHz. For the experiments, both instrument configurations employed the same free-space calibration procedure and controlled measurement setup. The result also compares well with the theoretical reflection coefficient. As expected, the response below 3 GHz does not merit consideration due to antenna gain limitations at lower frequencies. Additional reflection coefficient measurement of materials that have different dielectric properties was conducted with this free-space CMR-antenna prototype. A group of dielectric materials, including the silicone rubber sheet, a Skycast material (er = 12, thickness = 2.6cm), an FR4 board (er = 4.34, thickness = 0.2cm), and a Plexiglas material (er = 2.7, thickness = 2.4cm), was selected for measurements. The results, as summarized in Figure 3.8 and Figure 3.9, show good agreement with the theorectical values. 22 t 3 3CD <n c o Q. W CD DC 0.5 — i'n Ref short Offset short MUT Ac : fit : * ^ t l "yr^^lmll j i T T l >||ii<llil> 0 „ ( HIT CO £ o Q -0.5 0) 1 l': _: I i". : 11: i l | E Time-Gated & Calibrated - '-'• • i :: •l ~ i -1 - 0.5 1 Time (ns) 1.5 Figure 3.5 Time-gated and calibrated time-domain responses from the standards and MUT DQ -5 o f 1-15 o O c -20 o 1-25 1-30 -35. % 5 o CMR/antenna PNA/antenna Theory 6 7 8 9 Frequency (GHz) 10 11 12 Figure 3.6 Measured and calibrated magnitude of complex reflection coefficients of a silicone rubber sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMR-antenna system are compared with those obtained using commercial PNA-antenna setup as well as with the theoretical values. 23 o CMR/antennaPNA/antenna Theory 200 0 T3 ts c 100 0 o o O 0- ••G - 1 0 0 0 0 !^ -200 5 6 7 8 9 Frequency (GHz) 10 11 12 Figure 3.7 Measured and calibrated phase of complex reflection coefficients of a silicone rubber sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMR-antenna system are compared with those obtained using commercial PNA-antenna setup as well as with the theoretical values. 24 Fmynmcf |GH2| (a) I1 0 -10 i -20 J; f " |i * £ -30 '•} '* ii \ ;; ! ; '< ';f ~ :' i J i - -40 •SO ' ;:' f i —J i 1 — - i ~ - 7 B Ffetjfl?Has) (b) Figure 3.8 Measured reflection coefficients of (a) a Silicone rubber sheet (er = 3, thickness = 1.2cm) and (b) a Skycast material (er = 12, thickness = 2.6cm) using the free-space CMR-antenna prototype. Results (circles) compare favorably with the theory values (dotted lines). 25 Freq(GHz) (C) 6 7 8 Frequency (GMt) (d) Figure 3.9 Measured reflection coefficients of (c) an FR4 board (er = 4.34, thickness = 0.2cm) and (d) a Plexiglas material (er = 2.7, thickness = 2.4cm) using the free-space CMR-antenna prototype. Results (circles) compare favorably with the theoretical values (dotted lines). 26 Chapter 4 Concept of Acoustic and Microwave Hybrid Modality This chapter presents the concept of the acoustic and microwave hybrid modality for breast cancer detection. 4.1 Theory In the acoustic and microwave hybrid modality, an acoustic wave is applied to induce tissuedependent deformation in breast interior, and low-power microwave signals also illuminate the breast through an antenna array which in turn measures the scattered microwave signals. The induced tissue deformation modulates the scattered microwave signals, contributing to the microwave Doppler scattering response detected by the receiver. This modulation effect of the tissue deformation on microwave scattering is illustrated with a simple 1-D example. Assume that a halfspace dielectric interface, displacing at an acoustic frequency um, is illuminated by a continuous plane wave at a microwave frequency LO0, the scattered microwave field from the interface is then written as Es(z,t) = ^{Te-j20ld{t)e+j(}lZej^t} (4.1) where F is the reflection coefficient of the dielectric interface, /3i = y- is the microwave wavenumber in the medium, and the interface displacement d(t) is d(t) = Am sm(ujmt) (4.2) The acoustic-modulated scattered field is then expressed as +00 Es{t)\z=o = Tcos(u0t - 2/31Amsm(umt)) = F ^ . Jn{ . m Ai ) cos ({u0 + num)t) (4.3) 27 where Jn{x) is referred to as the Bessel function of the first kind of order n. Assume that the displacement is much smaller than the microwave wavelength, i.e. Am <C Ai which is true for the proposed application in this dissertation, and simplify the modulated scattered field using a small-argument Bessel approximation: t 1 JU(K) = { ±/c/2 for n =0 for \n\ = 1 0 for \n\ > 1 The acoustic-modulated scattered field is approximated as Ea(t)\z=0 2irA 2irA = Tcos(wot) + T — ^ cos((o;o + um)t) - T — ^ cos((u;o - uJm)t) Ai (4.4) Ai The first term in this expression is the fundamental scattered component which depends on the dielectric properties of the scatterer. The second and third terms represent the Doppler scattered components which depend on both dielectric and elastic properties of the scatterer. This theoretical analysis demonstrates the advantage of introducing acoustic modulation into microwave scattering, i.e., utilizing both dielectric and elastic properties contrasts in breast tissue instead of just one. 4.2 Methods of Inducing Tissue Deformation There are two main methods of coupling acoustic power into human tissue, as used in various biomedical applications: (a) a global excitation approach [50, 51, 52, 53, 54, 55], and (b) a local excitation approach [56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72]. For breast cancer detection, potential ways of inducing tissue deformation in breast are illustrated in Figure 4.1. External acoustic or mechanical excitation with a compression plate is a contact method. Local excitation with a focused acoustic transducer, however, can be non-contacting, distribute acoustic power deep in breast, and be configured for real-time beam steering as needed. Regarding the two options, the local excitation approach with the use of a focused acoustic transducer usually generates smaller displacement in human tissue (less than a fraction of mm), whereas the global excitation approach with the use of a compression plate can generate relatively larger deformation. The induced deformation in breast tissue directly affects the signal strength of the microwave 28 Doppler scattering. Numerical studies of the hybrid modality using both acoustic excitation approaches are presented in the following chapters. antenna compression plate _ f antenna breast j > - focused transducer (a) (b) Figure 4.1 Potential ways of inducing tissue deformation in breast: (a) global excitation with a compression plate, and (b) local excitation with a focused acoustic transducer 29 Chapter 5 Global-Acoustic and Microwave Hybrid Modality This chapter presents the numerical study of a global-acoustic and microwave hybrid modality which employs an external mechanical or acoustic excitation approach for inducing tissue deformation in breast. In addition, a synthetic means of tracking and removing bulk translation through the use of an antenna array and predictive signal processing is investigated. Finally, the limitations of using this global excitation approach in the hybrid modality are discussed. 5.1 2-D Numerical Simulation Modeling the acoustic (or mechanical) wave and microwave interaction in the hybrid modality requires a multiphysics numerical technique. Two distinct time scales are involved in this problem: the shorter time scale associated with high-frequency microwave illumination and the longer time scale associated with low-frequency acoustic (or mechanical) excitation. Simulating such wave interaction with a time-domain numerical technique such as the finite-difference timedomain (FDTD) method [73] would require a prohibitive number of time steps in computation. Thus, an efficient approach, the quasi-stationary (QS) approximation [74], is employed to solve this problem. In this multiphysics problem, the tissue target, under mechanical or acoustic excitation (~kHz), deforms sufficiently slow in comparison with the high-frequency (~GHz) microwave illumination. Thus, the microwave scattered fields should consist of, at any time t, the field pattern corresponding to the instantaneous deformed tissue target at any time Tj on the acoustic time scale, calculated as if the target is stationary. This field pattern is termed quasi-stationary. This QS method remains valid under the condition that ujm <C UJQ, where um is the angular frequency of the 30 low-frequency wave that induces tissue deformation and UJQ is the angular frequency of the highfrequency wave that interacts with the deforming tissue target which results in modulated scattered fields. Furthermore, the induced tissue deformation in breast is much smaller in comparison with the microwave wavelength. Resolving such small deformation in the numerical grid also results in a heavy computational load for FDTD. Thus, sheet boundary conditions (SBCs) [75, 76, 77] are adopted to efficiently model the microwave interaction with the deforming tissue target excited by acoustic (or mechanical) wave. With the SBC method, the excited object is modeled as a stationary object with time-varying sheet boundary conditions which account for the object perturbation on the relatively slow acoustic (or mechanical) time scale. SBCs were first postulated and derived for a 2-D dielectric object by Lawrence et al. [75], wherein SBCs were also implemented with the Method of Moments (MoM) in 2-D. More recently, Buerkle et al. [77] implemented SBCs with FDTD in 2-D. Overall, the numerical modeling of the global-acoustic and microwave hybrid modality involves two types of simulations. First, a mechanical simulation, based on a finite element method (FEM) solution of a plain-strain partial differential equation (PDE), is conducted to compute tissue deformation under mechanical excitation which is then imported into the subsequent electromagnetic simulation. Next, electromagnetic simulations, based on a FDTD solution of Maxwell's equations with the use of SBCs, model the microwave interactions with stationary snapshots of the the tissue target. In FDTD-SBC implementation, the total scattered field is decomposed into the unperturbed field scattered by the unperturbed object, E 0 , and a perturbation component which can be reconstructed as the Doppler scattering component, Ed, associated with the vibration of the object. The steps of the electromagnetic simulations are as follows: (1) An FDTD simulation is conducted to compute the scattered fields from the unperturbed object illuminated by the incident electromagnetic wave. (2) The equivalent surface currents for the SBCs are calculated using the unperturbed field solutions and peak boundary displacements computed from the mechanical simulation. Detailed expressions for the currents are given in [76] and are valid for small displacements 31 relative to the microwave wavelength — a regime that is indeed consistent with the proposed application with a wavelength of 1 cm in the tissue medium for the center microwave frequency 5 GHz. (3) FDTD-SBC simulations are conducted to compute perturbation components using equivalent surface currents on the stationary object as the only sources. (4) The perturbation components are then demodulated to reconstruct the magnitude of the Doppler scattering component. Implementation of the FDTD-SBC method for a 2D TM2 case is validated by comparing its results with those generated by a brute-force method in which the boundary of the scattering object is repositioned. This benchmark approach, referred as FDTD-REP, requires two FDTD simulations in which the scattered fields are computed for the same microwave illumination but different stationary scattering objects: an unperturbed object and the same object under peak deformation. The FDTD-REP solution, AEs(t), is calculated as the difference between the two sets of scattered field solutions. It can be shown analytically with the same 1-D problem described in Section 4.1 that AEs{t)\z=0 = rcos(o;ot) - Tcos(io0t = 1fixAm) r{cos(cu0t)(l — cos(—r—^)) — sin((j0i) sin(———)} Ai Ai (5.1) With the assumption that the displacement is much smaller than the microwave wavelength, i.e. A-m <C Ai, the expression can be approximated as 4-7T A AEa(t)\*=o = -r—-^sin(^t) (5.2) Ai Recall from Section 4.1 that the theoretically derived expression for the Doppler component of the acoustic-modulated microwave scattered field in this 1-D problem is Ed{t) = T cos((cu0 + ujm)t) - r Ai = —r——^ sin(u)0t) sm((dmt) Ai cos((o;o - wTO)t) M (5.3) 32 Then it is evident that the amplitude of AES is an excellent approximation to the amplitude of the Doppler component, provided that the displacement is much smaller than the microwave wavelength. Therefore, the use of the amplitude of the FDTD-REP solution as a benchmark is appropriate for the small deformation case considered here. The scattering object in the 2D TM2 validation study is a 4-mm-diameter cylinder (er = 50, as = 4.91 S/m at 5 GHz) in a homogeneous background material (er = 44, as = 4.16 S/m at 5 GHz). The radius of the inclusion is assumed to vary sinusoidally in time and uniformly across azimuthal angle, with a peak radial deformation of 0.45 mm. The incident wave is launched by an electric current line source 7.5 cm away from the inclusion center. Electric field observation points are 4 cm away from the inclusion center. Figure 5.1 shows the normalized amplitude of Ed and A E S computed at 5 GHz using FDTD-SBCs and FDTD-REP, respectively. The excellent agreement between these two set of results demonstrates the accuracy of the FDTD-SBC implementation. In the 2-D numerical testbed of the global-acoustic and microwave hybrid modality, a scenario that is challenging for a pure microwave detection scheme that is based solely on dielectric contrast is considered, i.e., investigating microwave scattering from a cylindrical inclusion of malignant tissue (er = 50, as = 4.91 S/m at 5 GHz) surrounded by normal tissue that has significant fibroglandular content (er = 35, as - 3.09 S/m at 5 GHz). For reference, scattering from a fibroglandular mass (er = 44, as = 4.16 S/m at 5 GHz) surrounded by the same normal tissue environment is also considered. In both cases, the scattering object (tumor or fibroglandular inclusion) is a 4-mm-diameter circular cylinder. The dielectric properties values assigned to the different tissues are derived from the data reported in [10]. While the tumor detection scenario described above involves a relatively small dielectric contrast between the scattering object and the surrounding medium, there is a large elastic properties contrast. The Young's modulus of the tumor is assigned as E = 12 kPa, while the surrounding normal tissue is modeled with the Young's modulus reported for fibroglandular tissue (E = 1.8 kPa) as reported in [30]. For the reference scenario, Young's modulus throughout the entire phantom (both the fibroglandular inclusion as well as the surrounding normal tissue) is assumed as E = 1.8 kPa. All tissue types are modeled using the same Poisson's ratio (y - 0.495). 33 -REP, [\Es|2/|Ei|2 SBC, |Ed|2/[Ei[2 -40 m T3 -50 -60 180 4> (deg) 360 Figure 5.1 Comparison of the amplitudes of the normalized Doppler component computed via FDTD-SBCs (dashed line) and the normalized difference between stationary scattered field solutions computed via FDTD-REP (solid line). The 2D scattering object is a deforming 4-mm-diameter circular cylinder undergoing a peak radial displacement of 0.45 mm. compression plate normal fibroglandular tissue 1 6L y" L=10cm — inclusion W=12cm Figure 5.2 Phantom configuration for the hybrid simulations, with axis of cylindrical inclusion into the page 34 In the mechanical simulation, the inclusion is embedded 5 cm below the upper surface of the 12 cmxlO cm domain, as shown in Figure 5.2. A mechanical compression plate is placed on the upper surface to apply an excitation with a peak strain of 5L. In simulation, applied peak strains from 1 mm to 10 mm (i.e., SL/L = 1-10%) have been tested. The structural mechanics module of Comsol Multiphysics software is used to conduct this FEM-based simulation and obtain the tissue-dependent deformation. The total deformation of the inclusion under mechanical compression consists of both elastic deformation and bulk translation. In order to isolate the tissue-dependent elastic deformation of the inclusion only, the bulk translation of the inclusion center 5L' is subtracted from the total simulated deformation. This calibrated deformation data from the mechanical simulation is imported into the FDTD-SBC electromagnetic simulation. The sampling density of the FDTD grid is N = Xi/Ax = 40, where Ai is the wavelength at 5 GHz in the surrounding normal tissue. The incident wave is launched by an electric current line source 7.5 cm away from the scatterer center. In this feasibility study, complicating effects due to scattering from the breast surface are avoided by filling the entire electromagnetic computational domain (aside from the inclusion region) with the normal breast tissue background medium. Electric field observation points are 4 cm away from the scatterer center. The mechanical simulation results, reported in Figure 5.3, show the elastic deformation of inclusion boundaries for the malignant tumor or normal fibroglandular inclusion, when a 10% peak strain applies. Bulk translation for both cases have been removed to separate out the tissue-dependent elastic deformation. The normal fibroglandular tissue inclusion experiences more than three times the boundary deformation of the malignant mass. The peak deformation of the top boundary of the inclusion (coordinates x' = 0, y' = 2 mm) is 0.24 mm for the normal fibroglandular inclusion and 0.07 mm for the malignant tumor. The subsequent electromagnetic simulation results, summarized in Figure 5.4, show the Doppler component in the backward ((f) = 180°) and forward (0 = 0°) scattered fields for the deforming scatterers as a function of the peak applied strain. Doppler scattering from the deforming normal fibroglandular inclusion is observed to be stronger than that from the deforming malignant tumor. The harmonic backscattered signal from the normal fibroglandular inclusion is 8 dB higher than \ 35 2. 1.5 - • (mm) 1 0.5 . \ ^»& * *' ^ * *^^ / V*j '•I 0 •if i il >% \ n -0.5 V* \» i original • *. 1 - ' - ' deformed (normal) - - - deformed (malignant) 11 Pi \A -1 "*4 -1.5 • ^^^^^^^J^^^^ -2 -2 -1 0 x ' (mm) • 1 2 Figure 5.3 Deformed inclusion boundaries, under a 10% peak applied strain (Bulk translation is removed). Peak deformation of boundary point (0, 2) is 0.24 mm for the normal fibroglandular inclusion and 0.07 mm for the malignant tumor inclusion. 36 -45 -55 CD 5. -65 • - -75Y ^ —©—normal, <> | = 180° X -85 V - • - normal, § = 0° * -95 0 1 malignant, i|) = 180° - * - malignant, $ = 0° 2 3 4 5 6 7 8 peak applied strain (%) 10 11 Figure 5.4 Computed Doppler components in the backward {<j> = 180°) and forward (</) - 0°) scattered fields under different peak applied strains, for malignant and normal fibroglandular scatterers. 37 that from the malignant tumor. In comparison, the backscattered signal from the unperturbed normal fibroglandular mass (in which case only the dielectric properties play a role) is 3.5 dB lower than that from the unperturbed malignant tumor. Hence, these results show the promise of using this hybrid modality to enhance the overall microwave scattering contrast for the purpose of distinguishing cancerous tumors from normal fibroglandular tissue in the breast. 5.2 Discussion of Limitations Despite showing some promising contrast enhancement results, the global excitation approach has some limitations. One assumption made in the 2-D numerical study is that the induced elastic deformation can be separated out from the total deformation by removing the bulk translation of the inclusion target. However, it is difficult for a practical antenna array to track the bulk translation of the target when external vibration applies to the breast. The bulk translation is not as tissue-dependent as the elastic deformation, and contributes to a majority of the total deformation. So with the bulk translation involved, there may not be much observed contrast in the microwave Doppler scattering between the malignant tumor and normal fibroglandular tissue cases. Therefore, in the following section, a synthetic means of tracking and removing the bulk translation through the use of an antenna array and predictive signal processing is investigated and tested with 3-D numerical simulations. 5.3 3-D Numerical Simulation with a Bulk Translation Removal Method As in 2-D, the QS approximation method is used to model this 3-D multiphysics problem. Two types of simulations are involved: (1) mechanical simulation and (2) electromagnetic simulation. However, for the preliminary 3-D work in this section, FDTD-REP method, as described in Section 5.1, is used to compute the amplitude of the Doppler scattering response from an object under mechanical excitation. In the 3-D numerical testbed as shown in Figure 5.5, a malignant tumor or a reference normal fibroglandular tissue clutter, both surrounded by normal breast tissue, is investigated. The same dielectric and elastic properties as those in the 2-D numerical testbed 38 in Section 5.1 are assigned to these phantom materials. The inclusion (tumor or fibroglandular tissue) is a 1-cm-diameter sphere. In the mechanical simulation testbed as shown in Figure 5.5(a), the inclusion is embedded 5 cm below the upper surface of the cylindrical phantom which has a height of 10 cm and a base diameter of 10 cm. A stainless steel compression plate is placed on the upper surface to apply an excitation with a 10% peak axial strain. The structural mechanics module of Comsol Multiphysics software is used to conduct this FEM-based simulation and obtain the tissue-dependent deformation. In the electromagnetic simulation testbed as shown in Figure 5.5(b), a multi-ring antenna array is modeled to transmit and receive microwave signals. The transmit antenna element is modeled as a A/2 dipole, and the receive antenna elements are modeled as electric field observation points. All the antenna elements are located 5 cm away from the central axis of the cylindrical phantom. The sampling density of the FDTD grid is N = Xi/Ax = 100, where Ai is the wavelength at 5 GHz in the surrounding tissue background. Similar to what is found in the mechanical simulation results in previous 2-D work, the total boundary deformation of the inclusion consists of both elastic deformation and bulk translation. The center of both the tumor and fibroglandular tissue targets experiences a 0.57 cm vertical bulk translation under the peak axial strain. This bulk translation is much larger than the elastic deformation of the inclusion which is on the order of a fraction of millimeter or less. The bulk translation is difficult to track and be removed by a physical antenna which would need to move synchronously with the breast phantom under the mechanical or acoustic excitation. Hence, a synthetic means of tracking and removing the bulk translation through the use of an antenna array and predictive signal processing is investigated as follow. Figure 5.6 and 5.7 illustrate simplified 2-D views of the 3-D phantom testbed, which helps explain the bulk translation removal procedure. First, as shown in Figure 5.6, stress Fi is applied to the upper surface of the phantom via the compression plate. A transmit antenna located at Zi launches microwave signals to the phantom, and the other antenna elements in the array receive the scattered microwave signals. Mechanical simulations are conducted to compute the tissue deformation under this applied stress. The spherical inclusion experiences both bulk translation and elastic deformation. Yet the elastic deformation is much smaller than the bulk translation. 39 stainless steel compression plate * -^ 5cm H=10em #- ¥ - D=10cm (a) <p=90° /C;;;;;;;;;;;;;;!!:* dipole 5. antenna ...•••' \ observation points (b) Figure 5.5 3-D numerical phantom testbed for the (a) mechanical simulation and (b) electromagnetic simulation of the global-acoustic and microwave hybrid modality. 40 Electromagnetic simulations then compute the microwave scattered fields from the inclusion both before and after the compression. We calculate the difference of the pre-compression and postcompression microwave scattered signals AE S for each receive antenna location. In this testbed, the tissue background is a homogeneous medium and the tissue target is a spherical inclusion. With these assumptions, at a certain receive antenna location Zn\, the EM wave propagation path from the transmit antenna Zt to the tissue target and then back to the receive antenna Zn\ is equivalently the same for the pre-compression and post-compression cases. Thus, the signals received by the antenna at Zn\ locations are approximately the same for pre-compression and post-compression, which results in a minimum value ('null') in the response AE S at the antenna location Zn\ of the array. Finding this minimum AE S response is referred as 'null' detection in the rest of the discussion. Next, we apply two different stresses, F\ and F2 as shown in Figure 5.7 to the tissue phantom, and calculate the bulk translation of the tissue target under compression based on the 'null' detection method as follow. The 'null' location in the AE S response is detected for each stress case. Then the distance L between the transmit antenna location Zi and the 'null' location Zn is obtained for each stress case. Assuming a linear stress-strain relationship of the elastic modulus of the tissue phantom, we then obtain the bulk translation Az of the tissue target under each applied stress as AZl = (L2 - Lx) • -^— (5.4) r2 — -f*i Az2^(L2-L1)--^-r (5.5) r2 — r\ After the bulk translation is approximated using the 'null detection' approach, the total deformation result of the inclusion from the mechanical simulation can be calibrated accordingly. The elastic deformation of the inclusion is obtained by subtracting the calculated bulk translation from the total deformation result. Figure 5.8 shows the elastic boundary deformation of the malignant tumor and the normal fibroglandular tissue clutter in this simulation testbed. The peak boundary deformation of the malignant tumor is 0.24 mm, in comparison with 0.73 mm of the normal fibroglandular tissue clutter case. The calibrated elastic deformation is then imported into the subsequent electromagnetic simulations. pre-compression post-compression stress F1 V +n* © i Az1 =-n1 -n1 Figure 5.6 An illustration of the 'null' detection procedure post-compression stress F1 post-compression stress F2 V T~* imAz1 L: V Az, =-n1 ^n2 gure 5.7 An illustration of the bulk translation calculation procedure 42 FDTD-REP simulations are conducted to compute the scattering solutions that approximate the magnitude of the Doppler component. Figure 5.9 shows the Doppler scattered fields received at different positions in the central antenna array which aligns with the inclusion center. <> / = 180° represents the dipole source location. A range of dielectric properties (e r =35~50, cr=4.162 S/m), representing the group of normal tissue most difficult to distinguish from tumor, is assigned to the normal fibroglandular inclusion in this testbed. Results show that the Doppler scattered field from the deforming malignant tumor is observed to be 5-10 dB weaker than that from the deforming normal fibroglandular tissue clutter. This Doppler scattering contrast is not only larger but also inverted in comparison with the fundamental scattering contrast of the pure microwave detection scheme. There is, however, a limit of the physical spacing between the antenna elements in a practical array. Errors in approximating the bulk translation using the 'null' detection method may be introduced with a coarse antenna spacing. Therefore, a linear estimation method is employed to interpolate and predict the microwave scattered signals received at locations between the physical antenna elements in the array. Consider an antenna array which is composed of N physical antenna elements evenly distributed in space. Let X{ (i=l,2,...N) represents the microwave scattered signal received by the ith physical antenna, and let X = [x1, x 2 , ...x N ], where the vector element in the matrix x1 = [xi,xi+1, ...XN,XI,X2, ...£J_I]T (i=l,2,...N). Also denote v = [vi,v2, . . . « M ] T as the microwave scattered signal at M locations between the physical antennas in the array. The goal of the linear estimation is to predict vt(i = 1,2, ...M) based on Xi(i = 1, 2, ...N) received by the physical antennas. The relationship of the detected and predicted signals can be represented as v = X« (5.6) where a is the coordinate vector of v with respect to X. Define R = X f f X and consider the inner product of v with X as j9 = X v (5.7) a = R~1f3 (5.8) The coordinate vector a can be obtained as 43 -,^-^-,—^p. unperturbed deformed normal deformed malignant ^HW.,.. -5 - 4 - 3 - 2 - 1 •..,jM«rr 0 1 x (mm) 2 3 Figure 5.8 Calibrated elastic boundary deformation of the malignant tumor and the normal fibroglandular tissue clutter. Bulk translation for both cases are removed with a synthetic tracking method through the use of an antenna array. The normal fibroglandular tissue clutter experiences about three times the boundary deformation of the malignant tumor. 44 -20 • malignant (e =50,0=4.910) normal (e=50,o=4.910) -30h m -a -40 Y normal (e : 35-50, a=4.162) LU **-«> - 5 0 LU -60 -70 u v%w ^*.-,^%fe ^ H 90 100 110 120 130 140 150 160 170 180 *(deg) Figure 5.9 The Doppler scattering response from the deforming malignant tumor and deforming normal fibroglandular tissues. Despite the dielectric properties variation in the normal fibroglandular tissue, the Doppler scattered field from the deformed malignant tumor is observed to be 5-10 dB weaker than that from the deformed normal fibroglandular tissue clutter. 45 -20 -30 m • error=0 •error=1 mm error=2 mm -40 normal fibroglandular clutter uT "50 10 W < en -60 malignant tumor -70h. -80 90 100 110 120 130 140 * (deg) 150 160 170 180 Figure 5.10 The Doppler scattering response from a deforming malignant tumor and a deforming normal fibroglandular tissue clutter, when different levels of errors in approximating the bulk translation are introduced in the testbed. 46 The microwave scattered fields at locations between the physical antennas can then be predicted based on the coordinate vector a. and signals received by the physical antennas X. Both the received and predicted scattered signals are used in the 'null' detection method to approximate the bulk translation of the inclusion. Different choices of physical antenna spacing are tested in simulation (Table 5.1). It is found that larger spacing between the physical antennas results in a larger error in the approximated bulk translation. It also shows that the spacing between the physical antennas needs to be smaller than 7 mm in order to obtain a relatively accurate bulk translation that has a moderate level of error, less than 2 mm. Errors introduced in the bulk translation approximation directly affects the computed Doppler scattering response in this computational testbed, as the total boundary deformation of the inclusion from the mechanical simulation is calibrated by subtracting the bulk translation before it is imported into the electromagnetic simulation. Figure 5.10 shows the change in the Doppler scattered fields from the malignant and normal fibroglandular inclusions when different levels of errors in approximating the bulk translation (error: 1 mm or 2 mm) are introduced. When such an error exists, the receive antennas do not track the inclusion center exactly and therefore the the signal strength of the received Doppler response decreases for both the malignant and normal cases. However, the Doppler scattering contrast between the malignant tumor and normal fibroglandular tissue clutter remains distinguishable when only moderate levels of such errors are present. 5.4 Discussion of Limitations As pointed out earlier in Section 5.1, the global-acoustic and microwave hybrid modality has its limitation due to the bulk translation issue, as it is difficult for a practical antenna array to track the bulk translation of the inclusion when external vibration applies to the breast. With the bulk translation involved, not much contrast in the microwave Doppler scattering between the malignant tumor and normal fibroglandular tissue clutter cases may be observed. In Section 5.3, a synthetic method of tracking and removing the bulk translation through the use of an antenna array and predictive signal processing is investigated and numerically tested with a 3-D phantom testbed. It is found that small spacing between the physical antennas in 47 Table 5.1 Introduced errors in calculating the bulk translation. antenna spacing A (mm) error in calculating the bulk translation (mm) 0.2 0 0.4 0 1 0 2 0 3 0.2 4 0.2 5 0.4 6 1.0 7 1.9 48 the array is required so that only a moderate level of error in approximating the bulk translation is introduced. However, this desired small spacing may cause other problems such as mutual coupling effects between the antennas in a practical setting. Furthermore, this bulk translation removal method does not apply for arbitrarily-shaped target cases, heavily relies on the homogeneity of the breast phantom, and may complicate the operating procedure of the imaging system based on this hybrid modality. Considering these combined limitations, a focused-acoustic and microwave hybrid modality is proposed which will be presented in the following chapters. 49 Chapter 6 Focused-Acoustic and Microwave Hybrid Modality This chapter introduces the focused-acoustic and microwave hybrid modality and provides an overview of the 3-D numerical techniques for simulating this multiphysics problem. 6.1 Local Excitation with a Focused Acoustic Transducer In a focused-acoustic and microwave modality, a focused acoustic transducer deposits energy to induce tissue deformation in a region of interest (ROI) in breast. Focused transducers have been employed in many biomedical applications, from low-power uses in sonography and various elastography technologies [60, 61, 62, 63, 64, 65, 66, 67] to high-intensity uses in hyperthermia [68], thermal ablation [69, 70], and drug delivery [71, 72]. The transducer can be either geometrically focused or electronically focused. The latter type is considered in this investigation, as it provides better depth of focus as well as flexibility in convenient beam-steering. 6.2 Overview of the 3-D Numerical Techniques The numerical modeling of the focused-acoustic and microwave hybrid modality involves two types of simulations as shown in Figure 6.1, which also employs a QS approximation method as introduced in Chapter 5. (1) 3-D acoustic-structural simulation A coupled acoustic-structural simulation, based on the FEM solutions of the acoustic propagation PDE (Equation 6.1), the stress-strain governing equations (Equation 6.2 and 6.3), and Newton's Law of motion (Equation 6.4), are conducted to compute the acoustically excited tissue 50 deformation. (In these equations, p is the acoustic pressure, p is the density of the medium, c is the longitudinal acoustic wave speed in the medium, UJ is the angular acoustic frequency, q and Q are the acoustic source terms, a and r are the normal and shear stresses respectively, e and 7 are the normal and shear strain respectively, E is the Young's modulus, v is the Poisson's ratio, u = [u, v, w]T is the deformation, and F is the applied force (N/m) [78].) This simulation is conducted using the acoustic-structure interaction module of the Comsol Multiphysics software. In the simulation, an acoustic wave is launched as the excitation source, and the excited tissue-dependent deformation is computed and exported to the subsequent electromagnetic simulation. v- (-- v p - q ^=Q <*x l - i / Gy V &z TXy Tyz T~xz E (1 + !/)(! --2v) V l-v v 0 0 0 ex v 0 0 0 e 0 0 0 tz 2 0 0 Ixy 0 l-2i/ 2 0 lyz 0 l-2i/ V V \-v 0 0 u 0 0 0 (6.1) pcz p 0 0 0 0 v du dx dv dy £z dw dz f-x e du Ixy 1 dv 2 v (6.2) Ixz (6.3) dy "^ dx dv lyz 1 dw dz ^ Txz 2 dy du 1 dw dz dx -pu \x = F (6.4) (2) 3-D electromagnetic simulation Electromagnetic simulations, based on a FDTD solution of Maxwell's equations (Equation 6.5 and 6.6) using SBCs, model the microwave interaction with the deforming tissue target and compute the Doppler scattered response. (In the equations, E and H are the electric and magnetic 51 fields, e, a, fj,, and a* are the dielectric properties of the medium, and J s and M s are the electric and magnetic source terms.) In this dissertation, the SBCs for a 3-D dielectric object are derived and SBC implementation with a 3-D FDTD algorithm is validated. ^! at = <9E ¥ _ I v x E - i ( M s + <7*H) (6.5) 1 „ 1,, = - V xH--(J l + aB) (6.6) (j, ^ Details of the two types of simulations are presented in the following two chapters respectively. Elastic and Acoustic Properties N Acoustic-Structural Mechanics Simulation Dielectric Properties Figure 6.1 Computational steps of simulating microwave scattering response from a deforming target under focused-acoustic excitation 52 Chapter 7 3-D Simulation of Focused-Acoustic Induced Tissue Deformation This chapter presents the coupled acoustic-structural simulation for computing the focusedacoustic induced tissue deformation. 7.1 Acoustic Beamformer Design and Modeling In the acoustic-structural simulation, the acoustic wave is launched by a focused acoustic transducer which is modeled based on an acoustic beamformer design presented as follow. An acoustic beamformer is a spatial filter that can be used with an array of acoustic transducer elements to focus energy at some desired location in a spatial field. Beamformer is comprised of one complex weight in each transducer element channel. The amplitude and phase in each channel is chosen to obtain constructive interference at the focus location and destructive interference elsewhere. The goal of the beamformer design is to maximize the energy deposited in a desired ROI in breast while minimizing energy deposited throughout the remainder of the breast region. Assume that the focused transducer is a square array transducer composed of N x N small circular pistons, the frequency response associated with the propagation from the nth piston element to the focus location is denoted as r by an(um, r), where um is the angular acoustic frequency in operation. Follow [79] and let the piston surface move uniformly with a speed of Uoelu'mt normal to the baffle, the analytical propagation response at a far-field location r(r, 6) is obtained as an(ujm,r) =j — U0kae 2Ji(kasm9) ka sin 6 (7.1) 53 where a is the radius of the piston, p is the density of the background medium, c is the velocity in the background medium, and k is the wavenumber. Also define a(r) = [ai(ujm,r),a2(uJm,r),...,aN2(iom,T)}T and let the TV2 x 1 vector h = [hi, h2,..., hN2]T be the concatenation of the weight vectors from each channel. When we apply the complex weights to the transducer elements, the observed signal at a location r is denoted as Y(r) as Y(r) = h H a ( r ) Then the energy received at the location r is qe(r) = h i / a ( r ) a / / ( r ) h Thus, the energy deposited in the ROI $ / is expressed as hHQfh, where Qf= / a(r)a*(r)dr The goal of the weights design is to maximize the total transmit energy deposited in the region $ f and can be written as h = arg max h hHh (7.2) h = maxeig{Qf} (7.3) One solution for Equation (7.2) is where the operator maxeig means selecting the eigenvector that corresponds to the maximum eigenvalue of the matrix Qf. Assume that the square transducer, composed of 21 x 21 circular piston elements, has an aperture size of 2 cmx 2 cm and is immersed in a homogeneous background medium (p = 1000 kg/m3, c = 1500 m/s) that mimics the normal fibroglandular breast tissue. The acoustic operating frequency is 150 kHz. The desired lateral focal plane is 2 cm away from the transducer aperture. Follow the design approach above, the transducer design is implemented using the acoustic-structure interaction module of Comsol Multiphysics Software. Figure 7.1 shows an achieved focused beam spot 54 in the (a) lateral and (b) axial focal planes. The focal size is roughly 1 cm x 1 cm in the lateral plane and 1 cm x 3.5 cm in the axial plane, which is a reasonable spot dimension for targeting a ROI to search for early stage tumors usually less than 2 cm. 7.2 Induced Deformation Results Then, with the use of this transducer as the excitation source, 3-D acoustic structural simulations are conducted to compute the induced tissue deformation excited. In the 3-D simulation testbed as shown in Figure 7.2, a 1-cm-diameter spherical tissue target surrounded by a homogeneous normal fibroglandular tissue background (E = 1.8 kPa, Poisson's ratio u = 0.495) is investigated. This target is either a malignant tumor (E = 12 kPa, v = 0.495) or a normal fibroglandular tissue clutter (E = 1.82 kPa, v = 0.495) for comparison. Coupled acoustic-structural simulations are conducted using the acoustic-structure interaction module of Comsol Multiphysics software, and reveal a roughly 3:1 contrast in the maximal radial deformation 5 between normal fibroglandular and malignant cases. Tumor deforms less as it is much stiffen Figure 7.3 and Figure 7.4 show the time-varying boundary deformation S of the malignant tumor and the normal fibroglandular tissue clutter respectively. Boundary deformation of the targets in both the lateral and axial focal planes are shown, ipiaterai is the angle with respect to the y-axis in the y-z lateral plane and ipaxiai is the angle with respect to the x-axis in the x-z axial plane. The group of solid-line curves represents the deformation for different snapshots on the acoustic time scale. It is observed that the largest boundary deformation occurs at ipaxiai = 0° for both targets due to the acoustic illumination pattern in ROI generated by the focused transducer. The largest boundary deformation is 0.07 mm for the malignant tumor case and 0.23 mm for the normal fibroglandular tissue clutter case. Furthermore, no bulk translation is observed from the resulted animations of the boundary deformation for both targets, which certainly shows an advantage of using the local instead of global excitation approach. It is critical to assess the safety of using this type of focused transducer to induce tissue deformation in human breasts. Applying high levels of acoustic power may cause cavitation-related bioeffects in breast. Therefore, investigation on the indicator for the potential for cavitation-related 55 -1.5 -1 -0.5 0 0.5 1 1.5 y (cm) Figure 7.1 Acoustic intensity (W/cm2) in the (a) lateral and (b) axial focal plane generated by a focused square transducer in a normal fibroglandular tissue background: The achieved focal region is roughly 1 cm x 1 cm x 3.5 cm. 56 1.5 cm acoustic force Rx \ antennas < i target 2 cm normal fibroglandular tissue background Tx / i / antenna 2 cm Figure 7.2 3-D computational testbed for simulating the microwave Doppler scattering response from a tissue target excited by a focused acoustic transducer. The target is either a malignant tumor or a normal fibroglandular tissue clutter, surrounded by normal breast tissues. 57 E E I to 45 90 135 180 Materal E E 90 135 180 V.axial Figure 7.3 Focused-acoustic induced time-varying radial deformation of a tumor boundary. ^lateral is the angle with respect to the y-axis in the y-z plane and ipaxtai is the angle with respect to the x-axis in the x-z plane. The group of solid-line curves represents deformation for different snapshots on the acoustic time scale. Its largest boundary deformation is 0.07 mm. 58 0.3 0.2 £ E 0.1 0 2 M -0.1 to -0.2 -0.3 0 45 90 135 180 ^lateral 45 90 135 180 ^axial Figure 7.4 Focused-acoustic induced time-varying radial deformations of a normal fibroglandular tissue clutter boundary, ipiaterai is the angle with respect to the y-axis in the y-z plane and ip axial is the angle with respect to the x-axis in the x-z plane. The group of solid-line curves represents deformation for different snapshots on the acoustic time scale. Its largest boundary deformation is 0.23 mm. 59 bioeffects, such as the mechanical index (MI) or thermal index (TI), in this focused-acoustic and microwave hybrid modality needs to be addressed in future work. 60 Chapter 8 3-D Simulation of Microwave Interaction with a Deforming Tissue Target This chapter presents the derivation of SBCs for a 3-D dielectric object and the validation of SBC implementation with a 3-D FDTD algorithm. 8.1 Derivation of 3-D SBCs SBCs were postulated for a penetrable object in [75] as h x (Ei - E 2 ) = n x h x (ZfHi + Z | H a ) (8.1) -n x (H x - H 2 ) = n x h x' (>?Ei + y 2 s E 2 ) (8.2) where Zf and Yf are the sheet impedances and admittances on the exterior (i = 1) or interior (i = 2) boundary of the object, n is the unit vector normal to the object boundary, E i and Hx are the electromagnetic fields outside the object, and E 2 and H 2 are the fields inside the object. If Zf = 0 and Yf = 0, the above expressions are simplified to the source-free (unperturbed) boundary conditions as h x (Eoi - E 0 2 ) = 0 (8.3) -h x (Hoi - H 0 2 ) = 0 (8.4) The sheet impedances and admittances for a 3-D dielectric object are derived by enforcing equivalence in the scattered field solutions for a perturbed (r = a + 8) dielectric sphere, and an 61 Figure 8.1 Sheet impedances Zs and admittances Ys for a 3-D dielectric object are derived by enforcing equivalence in the scattered field solutions for a) a perturbed (r = a + 8) dielectric sphere, and b) an unperturbed (r = a) dielectric sphere with SBCs applied at the boundaries. 62 unperturbed (r — a) dielectric sphere with SBCs applied at the boundaries, as illustrated in Figure 8.1. First, consider an unperturbed dielectric sphere (radius a) illuminated by a time-harmonic plane wave E 1 = xe~jkircose (8.5) H 1 = y^e-jkircose Z\ (8.6) where Z\ and k\ are the intrinsic impedance and wavenumber, respectively, of the background medium. As described in [80], the incident fields can be represented as a superposition of TM r 11 and TE r modes expressed in terms of potentials A$ and Frr* +oo A r = -^TaMkir)P^cos9) (8.7) + 0O Fj = ^J2anJn(k1r)P^(cos9) uZi , (8.8) n=l where u is the angular electromagnetic frequency, an = j n " ( „"^ , and Jn(x) = •y/^Jn+i(x). The eigenfunction solution for the scattered fields outside the sphere can be written as COS( Ars = —^V6„^(M^(cos0), r>a (8.9) F/ = ^Vcn^(fc1r)Pn1(cos0), r>a (8.10) n—l The eigenfunction scattering solution inside the sphere can be written in a similar manner as , +oo A , = S^YtdnJn(k2r)Pfcoa6), . r<a (8.11) r<a (8.12) . +oo s Fr = ^ V e „ J „ ( f c 2 r ) ^ ( c o s ^ , LUZ/2 —, n—l where Z2 andfc2are the intrinsic impedance and wavenumber, respectively, of the sphere. Then, the electric and magnetic field solutions can be determined from the above potentials using the following relations: 1 rj 2 63 = EQ Erh = 1 ld2Ar 1 1 dFr JLUfj,e r drdO e r sin 9 d(f> (8.14) _L 1 d2Ar | lldFr jujfxe r sin 9 drd<f>"*" e r 89 (8.15) Hr i jufie d2 dr 2 1 1 dAr (i r sin 9 d<fi #9 k2)Fr (8.16) 1 l<9 2 F r jw[ie r drd9 (8.17) 1 02FP — r sin 9 drd(f) [ir d9 + juj/ie 11 dA. Hj, = (8.18) Enforcing the continuity of the tangential components of the fields at r = a yields the following mode coefficients: -y/^2fHJn'(kia)Jn(k2a) , y/efiliH^ (kia)Jn(k2a) -yJe2HiJn(kia)Jn(k2a) y/eifj,2Jn(kia)Jn'{k2a) - ^ex\i2H\ (8.19) '{ha)Jn'{k2a) y/e2iJi\H\ '(kia)Jn'(k2a) ^e1fi2Jn'(k1a)Jn(k2a) 2 • ar, - ^/e^Hk £r( )'/ (kia)Jn{k2a) y/e2/j,iHk - ^jr2H\ r(2) dr.. = + + (8.20) (8.21) (ha)Jn(k2a) (ha) Jn'(k2a) (8.22) C-T). yJa-ilXiH^, '(ha)Jn'(k2a) - ^/e1fi2Hn (k1a)Jn(k2a) where ei,//i are the dielectric properties of the medium outside the sphere, and e2, fj,2 are the dielectric properties of the sphere. When the sphere is perturbed with a small radial amount 8, radius a in Equations (8.19) and (8.20) is replaced with a + S. Assuming that 5 is much smaller than the microwave wavelength and the sphere radius, we apply the first-order Taylor series expansion of each term about 5 = 0. The expressions then represent the perturbed mode coefficients bn and cn: -yftm. Jn'(ha) + hSJn"(ha) Jn(k2a) + +v/iT^i Jn(ha) + h8Jn'(ha) \ftm r(2)" Hk> (ha) + h5H(2) (ha) -y^iM2 m2)(ha) k25Jn'(k2a) - \Jn'(k2a) + k2SJn"(k2a) an \Jn(k2a) + 2} + h8Hkr(2)' (ha) k25Jn'(k2a) Jn'(k2a) + k25Jn"(k2a) (8.23) 64 V£2Mi Jn(ha) + h5Jn'(ha) +y/eifj,2 Jn'(ha) Jn'(k2a) + + hdJn"(ha) k25Jn"(k2a) Jn(k2a) + k26 Jn'(k2a) (8.24) Cn — v^I H{2)(kxa) + h5H{2)'(ha)~\ H{2)\ha) \jn'(k2a) + h5H{2)"(ha)} + k26Jn"(k2a) \jn(k2a) + k2SJn'(k2a) The SBCs (8.1) and (8.2) are also applied to the eigenfunction scattering solutions for the unperturbed dielectric sphere of radius a to derive alternative expressions for the coefficients in terms of Z{, Z2, Y{, and Y2. The scattered fields can be represented as a superposition of TM r and TE r modes. For the TM r mode, applying the sheet boundary conditions to the E$ and He field components at r = a yields the coefficient bn: -y/m [jn'(ha) - if±ZiJn(ha)] euT2 \jn(ha) - JB-YfJ^ha)] r(2) (ha) - ^Z{H(2\ha)\ y/eZjIiW* \Hk2)(ha) -^jT2 [jn(k2a) + - ^-Y{H{2) jw«i \jn'(k2a) \Jn(k2a) + (ha)\ JjfcYjJn'(k2a) + JUIC2 if±Z°2Jn(k2a) ^-Y2SJn'(k2a) \Jn'(k2a) + ^Z}jn(k2a) (8.25) Similarly, for the TE r mode, applying the sheet boundary conditions to the E^ and He field components at r = a yields the coefficient <v ~^m [Uha) +^m ^T, [Jn'(ha) - i^YfUha)} \H^(ha) ^JT2 - ^-ZtJn'(ha)] - ^-Z{H{2) r(2)' (ha) jw/n [jn'(k2a) + ^Y2°jn(k2a) [jn(M + \Jn'(k2a) + r(2) lit™ (ha) - ^Y{H%\ha) ^Z&'.foaj ™Y2°Jn(k2a) \Jn(k2a) + JUfl.2 -^-Z^(k2a (8.26) We equate the coefficient (8.23) with (8.25) to obtain the sheet impedances and admittances for the TM r mode, assuming that 5 is much smaller than a and simplifying kS^j (x) ks m2 i Jn(x) ^J^TT = kS [(^) 2 - l ] , where x = ka. In a similar way, we equate the coefficient (8.24) with (8.26) to obtain the sheet impedances and admittances for the TE r mode. Consider the superposition of the TM r and TE r modes. The sheet impedances and admittances are derived to be: Z{ = -jhZ.S (—? K ha> (8.27) 65 Z\ = +jk2Z25 Yl _ K k2aJ (8.29) z1 1 jk26 Yf (8.28) - <£>'' (8.30) where Z\ and Z2 are the intrinsic impedances of the medium outside the sphere and the sphere itself, respectively. We note that the first-order expressions of 8.27-8.30 exhibit a large degree of similarity to the first-order sheet impedance and admittance expressions derived in 2-D in [75, 76]. As discussed in [75, 76], the sheet impedances and admittances derived from the first-order SBCs (8.1) and (8.2) depend on the mode number n. Second-order sheet impedances and admittances, independent of the mode number n and thus applicable to an arbitrarily shaped object, are derived using the second-order SBCs presented below. h x (E x - E 2 ) = n x h x -n x (Hi — H 2 ) = n x n x "' H ' + ;§< s 'i Hi)+ '* H2+ ! (B2 i H2) 7iEi+ ( Ei)+72E2+ ( E2) ^ ^l l ^4 (8.31) (8.32) Jj is the surface tangent along the object boundary. We postulate that 77, B, 7, and A in 3-D maintain the same degree of similarity to the corresponding expressions in 2-D (reported in [75, 76]) as observed between the first-order sheet impedances and admittances in 3-D (equations 8.278.30) and 2-D (reported in [75, 76]). Thus, our postulated second-order expressions are as follows: -jhZiS (8.33) V2 = +Jk2Z2S (8.34) B' (8.35) ??i = Bs2 7is = ' fci jZ28 = ' k2 ' Zi 72 = + jk25 (8.36) (8.37) (8.38) 66 Ai = - (8.39) As — A2 — (8.40) k2Z2 These expressions account for the time-varying boundary deformation as 5 changes on the acoustic time scale. The second-order SBCs, independent of mode number n, can be applied to an arbitrarily shaped 3-D dielectric object with arbitrary perturbation. 8.2 3-D FDTD-SBC Implementation As described in Section 5.1, the total scattered field E from the deforming object is decomposed into the unperturbed field scattered by the unperturbed object Eo and a perturbation component E. For a perturbed, stationary object, E is at the same microwave frequency as the incident wave. For a vibrating object, E can be reconstructed as the microwave Doppler component Ed. The implementation of the SBCs with a 3-D FDTD algorithm is as follow. First, an FDTD simulation is conducted to compute the scattered fields from the unperturbed object illuminated by the incident electromagnetic wave. Then, the equivalent surface currents on the object boundary which result from the object perturbation are calculated as follows: J B = n x ( H 1 - H 3 ) = -n x n x 7iE 0 i + ^ ( A i — E o i ) +72E02 + ^ ( ^ 2 ^ E 0 2 ) —Ms = hx (Ei — E 2 ) = hxnx r?iH01 |(i^H 0 1 ) + ,2H02 + | ^ 2 | H 0 2 ) (8.41) (8.42) Let the unperturbed field solutions E 0 i = E 0 2 = E 0 and H 0 i = H 0 2 = H 0 , and insert the derived sheet impedances and admittances (8.33) to (8.40) into (8.41) and (8.42) to obtain d 1 n x nx ju(e1 - e2)6E0 - -os( ( Hi M s = hx fix juj(jii - M 2 ) £ H 0 - • 5 -(( OS Use the dielectric properties e, 1' -i d )-^-E0) H2 u os )—£^-H0) €\ 62 . W OS - e\ + -?1 and ^ = /^o for both the background medium (i = 1) and the deforming object (i = 2). The surface currents are then expressed as J s = ft x h x [joj(e[ - e'2)6E0 + (CTI - a2)6E0] 67 M , = —n x h x 1 ( eijw T + ff 1 4(4Ho) ds t'2ju+f'dsK Apply the Fourier Transform rules [81]. The time-domain expressions of the surface currents are n x n x « - 4)<^o = —n x n x . 1 -2ft e-! + (<n - <J2)8SQ 1 -2?t. e' 2 ' . 5 ,_ d _ , , , 'ds ds (8.43) (8.44) These ^-dependent surface currents account for the object perturbation. Next, FDTD-SBC simulations are conducted to compute the perturbation component E associated with each snapshot on the acoustic time scale, using the equivalent surface currents on the unperturbed object boundary as the only sources. These currents are introduced by adding the appropriate source terms in the FDTD update algorithm. The electric surface current or electric field in Equation (8.43) is composed of the x-, y-, and ^-directed components which are not collocated in a standard Cartesian FDTD grid. And for a curved object such as a sphere, taking the cross product on the electric field and the unit vector normal to the object boundary means that the f-directed current component is contributed by all the three x-, y-, and ^-directed electric fields. So is the y- and £-directed current components. Therefore, in order to obtain the <?0 value for calculating the electric current component at a surface location, all the neighboring electric fields residing on the same side of the boundary are averaged. For example, to obtain the £Q value for calculating a x-directed surface current that locates on the interior boundary, all the neighboring y- and zdirected electric fields which are also inside the object boundary are averaged. Then the electric current is calculated using the SQ value as in Equation (8.43). For calculating the magnetic surface current at a surface location, the derivative of the magnetic field with respect to the surface tangent in Equation (8.44) is calculated with a spatial central-difference method, using the difference of two neighboring magnetic field values devided by an approximated area. This approximated area is calculated using the radius, the latitude angle, and the azimuth angle of the neighboring field locations. The computed perturbation component changes as the surface current sources change on the acoustic time scale. 68 Finally, the computed perturbation component solution that varies on the acoustic time scale is demodulated to reconstruct the Doppler scattering component E j . This is a two-step demodulation process. The theoretical analysis of the Doppler reconstruction method is provided in Appendix A. The computed perturbation component from the FDTD-SBC simulation is a time-domain pulse response which corresponds to a wideband signal in the frequency domain, thus taking the Fourier transform of this perturbation component can extract its phasor component at the center microwave frequency associated with each time snapshot. This phasor component of the perturbation component which varies on the acoustic time scale is then demodulated via Fourier transform to obtain the Doppler scattering component. 8.3 Validation For validation purpose, a 3-D numerical testbed as shown in Figure 8.2 is considered. The scatterer in the validation study is a deforming malignant tumor (er -50,a = 4.91 at 5 GHz) in a normal fibroglandular tissue background (er = 35, a - 3.09 at 5 GHz). Multiple canonically deforming tumors of different shapes are chosen for the validation studies. The tumor is either perturbed and stationary or vibrating. The incident electromagnetic wave is launched by a A/2 dipole 2 cm away from the tumor center. 360 electric field observation points are selected to be 1.5 cm away from the tumor center in a 2-D plane. The sampling density of the FDTD grid is TV = Ai/Ax = 100, where Ai is the electromagnetic wavelength at 5 GHz in the normal fibroglandular tissue background, which results in a computational grid size of 141 x341 x441. For the perturbed, stationary tumor cases, the benchmark results are obtained by directly meshing the perturbed tumor and conducting a standard FDTD simulation. For the vibrating tumor cases, the vibrating tumor is modeled with two stationary snapshots on the acoustic time scale, i.e. the unperturbed tumor and the perturbed tumor under peak deformation. We then directly mesh the tumor for both snapshots and conduct standard FDTD simulations. With the assumption that the deformation is much smaller than the microwave wavelength for a simple canonically vibrating case, the difference of these two stationary solutions can be calculated and used as a good approximation of the magnitude of Doppler component. 69 normal fibroglandular tissue background _ = i. m ^ Rx antennas ',1.5 cm 2 cm Tx "antenna Figure 8.2 3-D numerical testbed for validating the FDTD-SBC implementation: simulate the microwave scattering response from a deforming malignant tumor surrounded by a normal fibroglandular tissue background -10 -15 •SBC • benchmark -20 •E, -25 -45 0 45 V (deg) 135 180 Figure 8.3 The computed total scattered fields from a slightly bloated spherical tumor with a uniform radial perturbation across its boundary, tp is the angle with respect to the y-axis in the y-z plane, tp = —90° is the monostatic angle. 70 DO M— T3 (D (D •+—» +-» ro o co "to o -50" • > -180 -135 - 9 0 > -45 • •0 45 ¥(deg) 90 135 180 Figure 8.4 The computed total scattered fields from a perturbed spherical tumor with a non-uniform radial perturbation across its boundary, tp is the angle with respect to the y-axis in the y-z plane, ip = —90° is the monostatic angle. -50' • ' ' • •-180 -135 -90 -45 ¥ (deg) 0 45 90 135 180 Figure 8.5 The computed total scattered fields from a perturbed ellipsoidal tumor with a non-uniform radial perturbation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, ij: = —90° is the monostatic angle. 71 Figure 8.3 shows the total scattered fields E 0 + E computed at 5 GHz from a slightly bloated spherical tumor with a uniform radial perturbation across its boundary. The total scattered fields in this and the following figures are all normalized by the scattered field received at the tumor center location but in a homogeneous normal fibroglandular tissue background (i.e., without the tumor existence). The size of the unperturbed tumor is 1 cm in diameter. The uniform radial perturbation is 5 = 0.3 mm. ip m this and the following figures is the angle with respect to the y-axis in the y-z plane, representing different field observation locations. The monostatic angle is ip = —90°. The FDTD-SBC result shows good agreement with the benchmark result. For reference, the total scattered field for the unperturbed sphere is also included to demonstrate that E 0 + E differs from EoFigure 8.4 shows the total scattered fields computed at 5 GHz from the same spherical tumor which is perturbed with a non-uniform radial perturbation across its boundary. The radial perturbation across the boundary is 5(p,4>,ip) = <5mcos(20)| cos^l, where 5m = 0.3 mm and (f> is the azimuth angle in the x-y plane. Figure 8.5 shows the total scattered fields computed at 5 GHz from an ellipsoidal tumor which is perturbed with a non-uniform radial perturbation across its boundary. 2 2 2 The boundary of this ellipsoidal tumor can be represented as ^ + j - + ^ < 1, where a0 = 5 mm, xl1 /2 z!2 bo = 4 mm, CQ = 3 mm. The perturbed boundary can be represented as-^- + ^- + -^- < 1, where a = 5.4 mm, b = 4.3 mm, c = 3.2 mm. There is good agreement between FDTD-SBC results and the benchmarks for both cases, which shows the validity of FDTD-SBC implementation for an object under arbitrary boundary perturbation. Figure 8.6 shows the computed Doppler scattered fields from a sinusoidally-vibrating spherical tumor with a uniform peak radial deformation across its boundary. The Doppler scattered fields in this and the following figures are all normalized by the scattered field received at the tumor center location but in a homogeneous normal fibroglandular tissue background (i.e., without the tumor existence). The size of the tumor is 1 cm in diameter. The uniform peak radial deformation is 5(t) = 5mcos(u>mt), where 5m = 0.3 mm and u>m is the angular acoustic frequency. For a simple canonically vibrating case such as this, FDTD-SBC simulation for only one snapshot on the acoustic time scale, i.e. when the object is under peak deformation, needs to be conducted to 72 compute the magnitude of Doppler component. Again, good agreement between the FDTD-SBC result and the benchmark is observed. Figure 8.7 shows the computed Doppler scattered fields from the same spherical tumor which is sinusoidally vibrating with a non-uniform peak radial deformation across its boundary. The non-uniform peak radial deformation is 5(p, cj), ip, t) = Sm cos(2<f>)\ costp\ cos(ujmt), where 5m 0.3 mm and oom is the angular acoustic frequency. The agreement between the FDTD-SBC and benchmark results in this case shows the validity of FDTD-SBC implementation for a vibrating object with an arbitrary deformation across its boundary. -?80 -135 -90 -45 0 45 ¥ (cleg) 90 135 180 Figure 8.6 The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor with a uniform peak radial deformation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, ip — —90° is the monostatic angle. 73 -30 CD U— T3 (D 0 a ro o w -40 SBC — benchmark • -50 \ -60 / V -70 i_ jy o. =80 o. o Q -90 -180 -135 -90 -45 0 45 ¥(deg) 90 135 180 Figure 8.7 The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor with a non-uniform peak radial deformation across its boundary, ip is the angle with respect to the y-axis in the y-z plane, ip = —90° is the monostatic angle. 74 Chapter 9 3-D Computational Investigation of Focused-Acoustic and Microwave Hybrid Modality Based on the numerical techniques presented in the previous two chapters, this chapter presents a 3-D computational investigation of the focused-acoustic and microwave hybrid modality. 9.1 Results - Microwave Scattering from a Canonically Deforming Target As a preliminary step, a canonically deforming target case in which the target is assumed to be deforming sinusoidally with a uniform peak deformation across its boundary is first investigated. Consistent with the acoustic-structural simulation testbed in Chapter 7, the tissue target under investigation is either a malignant tumor (er = 50, as = 4.91 S/m at 5 GHz) or normal fibroglandular tissue clutter. (er = 44, as = 4.16 S/m at 5 GHz). In both cases, the target (tumor or fibroglandular tissue clutter) is a 1-cm-diameter sphere surrounded by normal tissue that has significant fibroglandular content (er = 35, as = 3.09 S/m at 5 GHz). This represents a scenario that is challenging for the pure microwave detection scheme solely based on dielectric contrast to distinguish malignant tumor from normal tissue clutter. As shown in Figure 7.2, the incident wave is launched by a A/2 dipole 2 cm away from the target center. 360 electric field observation points are selected to be 1.5 cm away from the target center in a 2-D y-z plane. The sampling density of the FDTD grid is N = X\/Ax = 100, where Ai is the wavelength at 5 GHz in the normal fibroglandular tissue background, resulting in a computation grid size of 141 x341x441. A roughly 1:3 peak deformation contrast between malignant tumor and normal fibroglandular tissue clutter from the acoustic-structural simulations is observed as in the Chapter 7; thus, the tumor 75 boundary is assumed to be deforming sinusoidally with a uniform peak deformation of 0.1 mm, whereas the more elastic fibroglandular tissue clutter boundary is assumed to be deforming sinusoidally with a uniform peak deformation of 0.3 mm. For this simple canonically deforming case, FDTD-SBC simulation for only one snapshot on the acoustic time scale, i.e., when the target is under peak boundary deformation, is conducted to compute the magnitude of the Doppler scattering component. As reference, standard FDTD simulations are also conducted to compute the fundamental scattering fields from unperturbed targets (tumor or fibroglandular tissue clutter), as in the pure microwave detection approach. The simulation results are summarized in Figure 9.1 which shows the Doppler scattered fields from the two deforming malignant and normal tissue targets as well as the fundamental scattered fields from both unperturbed targets, lateral is the angle with respect to the y-axis in the y-z plane. Because the 1-cm-diameter spherical target under the microwave illumination in this testbed is a very weak back-scatterer, the strongest responses in both the fundamental and the Doppler scattering are observed to be at the forward scattered field location (ipiaterai = 90°). The magnitude of Doppler component from the deforming malignant tumor is observed to be lower than that from the deforming normal fibroglandular tissue clutter in the same normal tissue background. There is a 6-8 dB contrast in the Doppler signal strength between the two types of targets, which is observed to be larger and inverted in comparison with the contrast in their fundamental microwave scattering signals. 9.2 Results - Microwave Scattering from a Focused-Acoustic Induced Deforming Target Now the focused-acoustic excited tissue target case is addressed, wherein the simulated deformation results induced by the focused transducer, as shown in Figure 7.3 and Figure 7.4, is imported into the 3-D electromagnetic simulations. The orientation of the acoustic transducer and the microwave antennas is illustrated in Figure 7.2. The same dielectric properties as in the canonically deforming case are assigned to the targets investigated. FDTD-SBC simulations for a series of 12 snapshots on the acoustic time scale are conducted to compute the perturbation scattering component associated with each snapshot. The sampling number within one acoustic period is 76 -10 -20 •o » -30 a> "§ -40 i_ (U (U 8 -50 -60 -180 -135 -90 -45 0 45 de lateral < 9) 90 135 180 Figure 9.1 Simulated magnitude of Doppler scattering responses of a canonically deforming target, ipiaterai is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90° is the monostatic angle. There is a 6-8 dB contrast in the Doppler signal strength between the uniformly-deforming malignant tumor (5=0.1 mm) and normal fibroglandular tissue clutter (5=0.3 mm) cases. This contrast is larger and inverted, in comparison with the contrast in fundamental scattering from the two types of targets as in pure microwave detection techniques. 77 m T3 to O U) 0 -10 -20 -30 -40 -50 -60 -70 -80 -90 • normal fibroglandular clutter •malignant tumor -10' -180 -135 -90 fundamental st 1 Doppler -45 0 lateral 45 (cleg) 90 135 180 Figure 9.2 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced deforming target, ipiaterai is the angle with respect to the y-axis in the y-z plane and ^lateral = —90° is the monostatic angle. The Is* Doppler scattering contrast between the deforming malignant tumor (Smax=0.07 mm) and normal fibroglandular tissue clutter (5max=0.23 mm) cases is larger and inverted, in comparison with the contrast in the fundamental scattering from the two types of targets as in pure microwave detection techniques. 78 determined based on a trade-off of the requirement to reconstruct higher-order Doppler field components of interest and the computational economy. In our simulation, we are mainly interested in the 1st Doppler scattering component. Then we demodulate the perturbation phasor component that varies along the acoustic time scale via Fourier transform to reconstruct the Doppler scattering components. A theoretical analysis of the Doppler reconstruction method is provided in Appendix A. Figure 9.2 shows that the l 5 i Doppler scattering contrast between malignant tumor and normal fibroglandular tissue clutter is again larger and inverted, in comparison with the contrast in the fundamental scattering from the two types of targets as in pure microwave detection techniques, which offers potential for using this hybrid modality in enhancing the overall microwave scattering contrast for the purpose of distinguishing cancerous tumors from normal fibroglandular tissue clutter in breast. 9.3 Results - Impact of Dielectric and Elastic Properties Variations on the Microwave Doppler Scattering In the hybrid modality, the signal strength of the focused-acoustic modulated Doppler scattering depends on both dielectric and elastic properties in breast tissue. Therefore, with this 3-D focused-acoustic and microwave computational testbed, a series of simulations are conducted to investigate the impact of dielectric properties variation and elastic properties variation of the deforming target on the Doppler scattering responses. First, the dielectric properties of the target in this computational testbed are chosen to vary. The large-scale characterization data reported in [10] suggests that the dielectric properties of the normal fibroglandular tissue vary much more than those of the malignant tumor. Therefore, the following dielectric properties are assigned to the target under investigation: the 50th percentile data of normal fibroglandular tissue (er=44, CT=4.16 S/m), the 75th percentile data of normal fibroglandular tissue (er=48, cr=4.71 S/m), and the 50th percentile data of malignant tumor (er=50, cr=4.91 S/m). For all the three cases, the elastic properties of the target are assigned as (E=1.82 kPa, ^=0.495). All other parameters remain the same as in the previous focused acoustic and microwave simulation testbed. Results in Figure 9.3 show that the Doppler signal strength increases 79 with the increase of the dielectric properties of the target. Since the same elastic properties are assigned, the targets experience the same deforming behavior under the focused-acoustic excitation in all three cases. Therefore, the observed Doppler scattering contrast solely depends on the dielectric properties contrast of the targets. It is reasonable to conclude from this result that the Doppler scattering strength will become weaker for the group of normal breast tissue that has lower dielectric properties such as the adipose-dominated tissue. However, this type of tissue is much easier to distinguish from malignancy using the pure microwave detection scheme, as there is large dielectric contrast between the two tissue types. The main potential of the hybrid modality exists in situations such as distinguishing the normal breast tissue that has significant fibroglandular content from the malignant tumor. After examing variations of dielectric properties, the elastic properties of the target in this computational testbed are chosen to vary. The elastic properties of the malignant tumor vary widely depending on its cancerous type and stage, resulting in a range of elastic contrast between malignant tumor and normal fibroglandular tissue as reported in the literature [30, 31, 32, 33, 34]. Thus the following elastic properties are assigned to the target under investigation: E=12 kPa, 9 kPa, and 7.5 kPa, that represent the elastic properties of different types of malignant tumors, and E=l .82 kPa that represents the stiffness of a normal fibroglandular tissue clutter; v = 0.495 is assigned to all these target cases. The dielectric properties of the target under investigation is assigned as: e r =50, cr=4.91 S/m. All other parameters remain the same as in the previous focused-acoustic and microwave simulation testbed. In the acoustic-structural simulations, the largest radial deformation of the target across its boundary induced by the focused acoustic transducer is found to be: 0.07 mm (E=12 kPa target), 0.11 mm (E=9 kPa), 0.17 mm (E=7.5 kPa) and 0.23 mm (E=1.82 kPa). The deformation of each target case is then imported to the subsequent electromagnetic simulations to compute the Doppler scattering response. The targets, assigned with different elastic properties, experience different deforming patterns under the focused-acoustic excitation, which results in different scattering patterns across the observation angle as shown in Figure 9.4. However, the averaged Doppler signal strength detected across the observation angles is found to be: -78.6 dB (E=12 kPa target), -75.8 dB (E=9 kPa), -71.9 dB (E=7.5 kPa), and -66.4 dB (E=1.82 kPa), which 80 shows that the overall Doppler signal strength also increases with the decrease of the elastic modulus of the target. The result also implies that this hybrid modality performs better in distinguishing stiffer breast carcinomas from normal fibroglandular tissue clutter. Overall, the computational investigation of the hybrid modality demonstrates a larger and inverted microwave Doppler scattering contrast between the malignant and normal fibroglandular tissues in comparison with their fundamental scattering contrast of a pure microwave detection scheme. This reveals the potential of the hybrid modality for improved sensitivity in detecting a malignant tumor in a heterogeneously dense breast environment. 81 m T3 2 CD Q. O Q 180 Figure 9.3 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced deforming target which varies in its dielectric properties. The assigned dielectric properties are: the 50th percentile data of normal fibroglandular tissue (er=44, cr=4.16), the 75th percentile data of normal fibroglandular tissue (er=48, CT=4.71), and the 50th percentile data of malignant tumor (er=50, (7=4.91). ipiaterai is the angle with respect to the y-axis in the y-z plane and ^lateral = —90° is the monostatic angle. 82 -30 CQ T> -nu <> / -50 •a <u ip - - - e =50,a=4.91 r e =50,c=4.91 r e =50,a=4.91 r e=50,o=4.91 r • • T) P CI) -KO S/m, E=12 kPa S/m, E=9 kPa S/m, E=7.5 kPa S/m, E = 1 . 8 2 k P a (0 u -70 ' ^ i_ \ CD n a. o Q -HO • r S ^ - , ' ' — 1 «. /,'' -. jt~ .*'.*..'*'— '+ + ' v " - ^ S _ _ ^ ^ ^ ^—"^ s - s -90 -100 -180 -135 ,• 'v.'- .-V -90 -45 .-''><-"-*».' '''•., .•• y s ••• -....•' /• .' Vs .-;' ' >~ \ V ~--' 45 , (deg) 90 135 180 Figure 9.4 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced deforming target which varies in its elastic properties. The assigned elastic properties are: E=12 kPa, 9 kPa, and 7.5 kPa that represents different types of malignant tumors, and E=1.82 kPa that represents a normal fibroglandular tissue clutter, v — 0.495 is assigned to all these target cases. ^lateral is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90° is the monostatic angle. 83 Chapter 10 Conclusion and Future Work This dissertation addresses two important issues in current microwave imaging systems for breast cancer detection: cost reduction and sensitivity improvement. The dissertation first presents the development of a low-cost CMR-antenna subsystem. The specific contributions include: (1) Cascaded mixer-incorporated bridge-T attenuators are used as the wave separators in this refiectometer design. This approach allows significant cost and size reduction of a microwave refiectometer. (2) A three-step free-space calibration method is developed for the subsystem. Unlike the conventional mechanical calibration procedure, this contactless method eliminates the need for connecting mechanical standards repeatedly at the reference planes and compensates for the antenna reverberation, thus providing better measurement repeatability. The method also allows flexibility in the choice of reference plane positions. (3) The free-space calibrated measurement results using the CMR-antenna prototype compare favorably with those using a VNA-antenna setup. Therefore, multiple sets of the reflectometerantenna subsystem can be used to construct a low-cost array-based hardware system for microwave breast imaging. Following this track, future research involves the development of a fully-populated receiver with the use of this free-space calibrated CMR-antenna subsystem. Specific tasks of potential interest include the following: 84 (1) Conduct simulation to investigate the mutual coupling effects between neighboring antenna elements in the array. Better understanding of the mutual coupling effects aids in the design of a fully-populated configuration. (2) Develop free-space calibration standards for use in a specific array. Depending on the array configuration (circular, planar, etc.), different free-space calibration standards will be needed. (3) Construct a fully-populated CMR-antenna array prototype, and investigate its performance in measuring scattered response in both time and frequency domains. This dissertation also presents an acoustic and microwave hybrid modality that exploits both dielectric and elastic properties contrasts in breast tissue. The specific contributions include: (1) Two methods of inducing tissue deformation in breast, global excitation with a compression plate and local excitation with a focused acoustic transducer, are investigated. (2) 3-D numerical techniques are developed to efficiently model this multi-physics problem. The FEM-based acoustic-structural simulations are conducted to compute breast tissue deformation excited either by mechanical vibration or by focused-acoustic illumination. Simulation results show that the malignant tumor which has higher elastic modulus deforms less than the normal fibroglandular tissue target. The FDTD based numerical method with the use of sheet boundary conditions is adopted and extended to 3-D to simulate the microwave scattering response of the acoustically excited tissue. (3) Despite showing some promising initial results, the numerical study of the global-acoustic and microwave modality reveals the limitations of employing an external excitation approach, including the difficulty of tracking and removing the bulk translation of the excited target via a physical or synthetic antenna array. (4) The numerical study of the focused-acoustic and microwave modality, however, shows that it avoids the bulk translation problem. The computational investigation demonstrates a larger and inverted microwave Doppler scattering contrast between the malignant and normal fibroglandular tissues in comparison with their fundamental scattering contrast of a pure microwave detection scheme. While the pure microwave imaging scheme creates a breast image of adipose-dominated and fibroglandular tissues, the hybrid modality can provide additional contrast in the microwave 85 Doppler scattering response to help distinguish the malignant from normal fibroglandular tissues and provide important diagnostic information. Therefore, the main role of this hybrid modality is in detecting malignant tumors in a heterogeneously dense breast. Following this track, future research tasks are as follows: (1) Simulate the focused-acoustic modulated microwave scattering response with a realistic 3-D breast phantom testbed, wherein a heterogeneous breast interior (preferably an anatomically realistic breast interior), skin layer, and coupling media layers for the acoustic transducer are included in the simulation testbed. (2) Perform further parametric studies of the focused-acoustic modulated microwave Doppler scattering strength. Parameters of interest include the size, the shape, and the density of the lesion, the properties of the heterogeneous tissue background, the operating frequency and power of the acoustic source, and the operating microwave frequency. (3) Conduct safety assessment of the focused-acoustic excitation approach employed in the hybrid modality. 86 LIST OF REFERENCES [1] P. T. Huynh, A. M. Jarolimek, and S. Daye, "The false-negative mammogram," Radiographics, vol. 18, no. 5, pp. 1137-1154, 1998. [2] V. P. Jackson, R. E. Hendrick, S. A. Reig, and D. B. Kopans, "Imaging of the radiographically dense breast," Radiology, vol. 188, pp. 297-301, 1993. [3] J. G. Elmore, M. B. Barton, V. M. Moceri, S. Polk, P. J. Arena, and S. W. Fletcher, "Ten-year risk of false positive screening mammograms and clinical breast examinations," New Engl. J. Med., vol. 338, no. 16, pp. 1089-1096, 1998. [4] M. Sabel and H. Aichinger, "Recent developments in breast imaging," Phys. Med. Biol, vol. 41, pp. 315-368, 1996. [5] S. H. Heywang-Kobrunner, "Nonmammographic breast imaging techniques," Current Opinion in Radiology, vol. 4, pp. 146-154, 1992. [6] A. J. Surowiec, S. S. Stuchly, J. R. Barr, and A. Swamp, "Dielectric properties of breast carcinoma and the surrounding tissues," IEEE Trans. Biomed. Eng., vol. 35, pp. 257-263, April 1988. [7] W. T. Joines, Y. Zhang, C. Li, and R. L. Jirtle, "The measured electrical properties of normal and malignant human tissues from 50 to 900 MHz," Med. Phys., vol. 21, pp. 547-550, April 1994. [8] S. S. Chaudhary, R. K. Mishra, A. Swarup, and J. M. Thomas, "Dielectric properties of normal and malignant human breast tissues at radiowave and microwave frequencies," Indian J. Biochem. Biophys., vol. 21, pp. 76-79, February 1984. [9] A. M. Campbell and D. V. Land, "Dielectric properties of female human breast tissue measured in vitro at 3.2 GHz," Phys. Med. Biol, vol. 37, pp. 193-210, 1992. [10] M. Lazebnik, D. Popovic, L. McCartney, C. B. Watkins, M. J. Lindstrom, J. Harter, S. Sewall, T. Ogilvie, A. Magliocco, T. M. Breslin, W. Temple, D. Mew, J. H. Booske, M. Okoniewski, and S. C. Hagess, "A large-scale study of the ultrawideband microwave dielectric properties of normal, benign and malignant breast tissues obtained from cancer surgeries," Phys. Med. Biol, vol. 52, pp. 6093-6115, 2007. 87 [11] S. Y. Semenov, A. E. Bulyshev, A. E. Souvorov, R. H. Svenson, Y. E. Sizov, V. Y. Borisov, V. G. Posukh, I. M. Kozlov, A. G. Nazarov, and G. P. Tatsis, "Microwave tomography:Theoretical and experimental investigation of the iteration reconstruction algorithm," IEEE Trans. Microw. Theory Tech., vol. 46, February 1998. [12] A. E. Souvorov, A. E. Bulyshev, S. Y Semenonv, R. H. Svenson, A. G. Nazarov, Y E. Sizov, and G. P. Tatsis, "Microwave tomography: A two-dimensional Newton iterative scheme," IEEE Trans. Microw. Theory Tech., vol. 46, November 1998. [13] P. M. Meaney, M. W. Fanning, D. Li, S. P. Poplack, and K. D. Paulsen, "A clinical prototype for active microwave imaging of the breast," IEEE Trans. Microw. Theory Tech., vol. 48, pp. 1841-1853, November 2000. [14] J. D. Zaeytijd, A. Franchois, C. Eyraud, and J. M. Geffrin, "Full-wave three-dimensional microwave imaging with a regularized Gauss-Newton method - theory and experiment," IEEE Trans. Antennas Propag., vol. 55, no. 11, pp. 3279-3292, 2007. [15] P. Kosmas, J. D. Shea, B. D. Van Veen, and S. C. Hagness, "Three-dimensional microwave imaging of realistic breast phantom via an inexact Gauss-Newton algorithm," in Proc. Antennas Propag. Society Intl. Symp., July 2008. [16] D. W. Winters, J. D. Shea, P. Kosmas, B. D. Van Veen, and S. C. Hagness, "Three-dimensional microwave breast imaging: Dispersive dielectric properties estimation using patient-specific basis functions," IEEE Trans. Med. Imag., vol. 28, pp. 969-981, July 2009. [17] X. Li and S. C. Hagness, "A confocal microwave imaging algorithm for breast cancer detection," IEEE Microw. Wireless comp. Lett., vol. 11, March 2001. [18] 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 Trans. Biomed. Eng., vol. 49, no. 8, pp. 812-822, 2002. [19] S. K. Davis, E. J. Bond, X. Li, S. C. Hagness, and B. D. Van Veen, "Microwave imaging via space-time beamforming for early detection of breast cancer: Beamformer design in the frequency domain," 7. of Electromagnetic Wave and Applications, vol. 17, no. 2, 2003. [20] S. K. Davis, H. Tandradinata, S. C. Hagness, and B. D. Van Veen, "Ultrawideband microwave breast cancer detection: A detection-theoretic approach using the generalized likelihood ratio test," IEEE Trans. Biomed. Eng., vol. 52, no. 7, 2005. [21] X. Li, S. K. Davis, S. C. Hagness, D. W. van der Weide, and B. D. Van Veen, "Microwave imaging via space-time beamforming: experimental investigation of tumor detection in multilayer breast phantom," IEEE Trans. Microw. Theory Tech., vol. 51, pp. 1856-1865, 2004. 88 [22] J. M. Sill and E. C. Fear, "Tissue sensing adaptive radar for breast cancer detection - experimental investigation of simple tumor models," IEEE Trans. Microw. Theory Tech., vol. 53, November 2005. [23] A. Wangsanata, "A planar vector network analyzer," Master's thesis, University of WisconsinMadison, 2002. [24] G. F. Engen, "The six-port reflectometer: An alternative network analyzer," IEEE Trans. Microw. Theory Tech., vol. 25, pp. 1075-1080, 1977. [25] G. F. Engen, "A (historical) review of the six-port measurement technique," IEEE Trans. Microw. Theory Tech., vol. 45, pp. 2414-2417, 1997. [26] F. R. de Sousa, B. Huyart, and R. C. S. Freire, "Low cost network analyzer using a six-port reflectometer," in Proc. IEEE MTT-S Intl. Microw. Optoelectronics Conf., 2001. [27] F. Wiedmann, B. Huyart, E. Bergeault, and L. Jallet, "New structure for a six-port reflectometer in monolithic microwave integrated-circuit technology," IEEE Trans. Instrum. Meas., vol. 46, pp. 527-530, 1997. [28] B. Huyart, F. Wiedmann, L. Jallet, E. Bergeault, R. Benelbar, and R. G. Bosisio, "Microwave measurement using Wheatsone's bridges," in Proc. European Microw. Conf., 1994. [29] R. A. Marsland, C. J. Madden, D. W. van der Weide, M. S. Shakouri, and D. M. Bloom, "Monolithic integrated circuits for mm-wave instrumentation," in Proc. IEEE Annual GaAs ICSymp., 1990. [30] A. Samani, J. Bishop, C. Luginbuhl, and D. B. Plewes, "Measuring the elastic modulus of ex vivo small tissue samples," Phys. Med. Biol., vol. 48, pp. 2183-2198, 2003. [31] A. Samani, J. Zubovits, and D. Plewes, "Elastic moduli of normal and pathological human breast tissues: An inversion-technique-based investigation of 169 samples," Phys. Med. Biol, vol. 52, pp. 1565-1576, January 2007. [32] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, "Viscoelastic shear properties of in vivo breast lesions measured by MR elastography," Magnetic Resonance Imaging, vol. 23, pp. 159-165, 2005. [33] A. R Sarvazyan, "Shear acoustic properties of soft biological tissues in medical diagnostics," /. Acoust. Soc. Am. Proc, vol. 93, pp. 2329-2330, April 1993. [34] T. A. Krouskop, T. M. Wheeler, F. Kallel, B. S. Garra, and T. Hall, "Elastic moduli of breast and prostate tissue under compression," Ultrasound Imaging, vol. 20, pp. 260-274, 1998. 89 [35] A. Thomas, S. Kummel, F. Fritzsche, M. Warm, B. Ebert, B. Hamm, and T. Fischer, "Realtime sonoelastography performed in addition to B-mode ultrasound and mammography: Improved differention of breast lesions?," Acad. Radiology, vol. 13, pp. 1496-1504, December 2006. [36] E. S. Burnside, T. J. Hall, A. M. Sommer, G. K. Hesley, G. A. Sisney, W. E. Svensson, J. P. Fine, J. Jiang, and N. J. Hangiandreou, "Differentiating benign from malignant solid breast masses with US strain imaging," Radiology, vol. 245, no. 2, pp. 401^10, 2007. [37] R. Sinkus, J. Lorenzen, D. Schrader, M. Lorenzen, M. Dargatz, and D. Holz, "High-resolution tensor MR elastography for breast tumor detection," Phys. Med. Biol, vol. 45, pp. 1649-1664, 2000. [38] A. Alizad, M. Fatemi, D. H. Whaley, and J. F. Greenleaf, "Application of vibro-acoustography for detection of calcified arteries in breast tissue," J. Ultrasound Med., vol. 23, pp. 267-273, 2004. [39] R. A. Kruger, K. K. Kopecky, A. M. Aisen, D. R. Reinecke, G. A. Kruger, and W. L. Kiser, "Thermoacoustic CT with radio waves: A medical imaging paradigm," Radiology, vol. 211, pp. 275-278, 1999. [40] R. A. Kruger, K. Stantz, and W. L. Kiser, "Thermoacoustic CT of the breast," in Proc. SPIE, vol. 4682, pp. 521-525, 2002. [41] G. Ku, B. D. Fomage, X. Jin, M. Xu, K. K. Hunt, and L. V. Wang, "Thermoacoustic and photoacoustic tomography of thick biological tissues toward breast imaging," Technol. Cancer Res. Treatment, vol. 4, pp. 559-565, 2005. [42] M. Xu, K. Geng, J. Xing, L. V. Wang, B. D. Fomage, and K. K. Hunt, "Breast cancer imaging by microwave-induced thermoacoustic tomography," in Proc. SPIE, vol. 5697, pp. 45-48, 2005. [43] D. M. Pozar, Microwave Engineering. Readings, Mass.: Addison-Wesley, 1990. [44] D. K. Ghodgaonkar, V. V. Varadan, and V. K. Varadan, "A free-space method for measurement of dielectric constants and loss tangents at microwave frequencies," IEEE Trans. Instrum. Meas., vol. 37, pp. 789-793, June 1989. [45] D. K. Ghodgaonkar, V. V. Varadan, and V. K. Varadan, "Free-space measurement of complex permittivity and complex permeability of magnetic materials at microwave frequencies," IEEE Trans. Instrum. Meas., vol. 39, pp. 387-394, April 1990. [46] N. Gagnon, J. Shaker, L. Roy, A. Petosa, and P. Berini, "Low-cost free-space measurement of dielectric constant at Ka band," in Proc. IEE Microw. Antennas Propag., vol. 151, pp. 271— 276, June 2004. 90 [47] M. Nakhkash, Y. Huang, W. Al-Nuaimy, and M. Fang, "An improved calibration technique for free-space measurement of complex permittivity," IEEE Trans. Geo. Remote Sens., vol. 39, pp. 453-455, February 2001. [48] F. C. Smith, B. Chambers, and J. C. Bennett, "Calibration techniques for free space reflection coefficient measurements," in Proc. IEE A-Science Meas. Tech., vol. 139, pp. 247-253, September 1992. [49] X. Li, S. C. Hagness, M. K. Choi, and D. W. van der Weide, "Numerical and experimental investigation of an ultrawideband ridged pyramidal horn antenna with curved launching plane for pulse radiation," IEEE Antennas Wireless Propag. Lett, vol. 2, pp. 259-262, 2003. [50] K. J. Parker, S. R. Huang, R. A. Musulin, and R. M. Lerner, "Tissue response to mechanical vibration for sonoelasticity imaging," Ultrasound Med. Biol, vol. 16, pp. 241-246, 1990. [51] Y. Yamakoshi, J. S ato, and T. Sato, "Ultrasonic imaging of the internal vibration of soft tissue under force vibration," Ultrasound Ferroelectr. Freq. Control, vol. 37, pp. 45-53, 1990. [52] R. Muthupillai, D. J. Lomas, P. J. Rossman, J. F. Greenleaf, A. Manduca, and R. L. Ehman, "Magnetic resonance elastography by direct visulization of propagating acoustic strain waves," Science, vol. 269, pp. 1854-1857, 1995. [53] S. Catheline, F. Wu, and M. Fink, "A solution to diffraction biases in sonoelasticity: The acoustic impulse technique," /. Acoust. Soc. Am., vol. 105, pp. 2941-2950, 1999. [54] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, "Shear modulus imaging with 2-D transient elastography," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 49, pp. 426^435, 2002. [55] M. Zhao, J. D. Shea, S. C. Hagness, D. W. van der Weide, B. D. Van Veen, and T. Varghese, "Numerical study of microwave scattering in breast tissue via coupled dielectric and elastic contrasts," IEEE Antennas Wireless Propag. Lett., vol. 7, pp. 247-250, 2008. [56] K. R. Nightingale, P. J. Kornguth, and G. E. Trahey, "The use of acoustic streaming in breast lesion diagnosis: A clinical study," Ultrasound Med. Biol., vol. 25, pp. 75-87, September 1998. [57] S. Calle, J. P. Remenieras, O. B. Matar, M. E. Hachemi, and F. Patat, "Temporal analysis of tissue displacement induced by a transient ultrasound radiation force," J. Acoust. Soc. Am., vol. 118, no. 5, 2005. [58] T. Sugimoto, S. Ueha, and K. Itoh, "Tissu hardness measurement using the radiation force of focused ultrasound," in Proc. IEEE Intl. Ultrasounics Symp., pp. 1377-1380, 1990. [59] W. F. Walker, "Internal deformation of a uniform elastic solid by acoustic radiation force," J. Acoust. Soc. Am., vol. 105, pp. 2508-2518, 1999. 91 [60] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov, "Shear wave elasticity imaging: A new ultrasonic technology of medical diagnostics," Ultrasound Med. Biol, vol. 24, no. 9, pp. 1419-1435, 1998. [61] E. A. Barannik, "Acoustic remote palpation: The influence of tissue viscosity on the excitation and relaxation of local impulse shear strain," in Proc. WCU, pp. 7-10, September 2003. [62] K. Nightingale, M. S. Soo, R. Nightingale, and G. Trahey, "Acoustic radiation force impulse imaging: in vivo demonstration of clinical feasibility," Ultrasound Med. Biol, vol. 28, no. 2, pp. 227-235, 2002. [63] J. Bercoff, M. Tanter, S. Chaffai, and M. Fink, "Ultrafast imaging of beamformed shear wave induced by the acoustic radiation force," in Proc. IEEE Int. Ultrasonics Symp., 2002. [64] E. E. Konofagou, M. Ottensmeyer, S. Agabian, S. L. Dawson, and K. Hynynen, "Estimating localized oscillatory tissue motion for assessment of the underlying mechanical modulus," Ultrasonics, vol. 42, pp. 951-956, 2004. [65] M. Fatemi and J. F. Greenleaf, "Ultrasound-stimulated vibro-acoustic spectrography," Science, vol. 280, pp. 82-85, 1998. [66] S. Calle, J. P. Remenieras, O. B. Matar, M. Defontaine, and F. Patat, "Application of nonlinear phenomena induced by focused ultrasound to bone imaging," Ultrasound Med. Biol, vol. 29, pp. 465^72, 2003. [67] R. Souchon, R. Salomir, O. Beuf, L. Milot, D. Grenier, D. Lyonnet, J. Y. Chapelon, and O. Rouviere, "Transient MR elastography (t-MRE) using ultrasound radiation force: Theory, safety, and initial experiments in vitro," Magn. Res. Med., vol. 60, pp. 871-881, October 2008. [68] C. J. Diederich and K. Hynynen, "Ultrasound technology for hyperthermia," Ultrasound Med. Biol, vol. 25, pp. 871-887, July 1999. [69] E. Kennedy, G. R. T. Haar, and D. Cranston, "High intensity focused ultrasound: Surgery for the future?," Br. J. Radiol., vol. 76, pp. 590-599, September 2003. [70] G. ter Haar, "Ultrasound focal beam surgery," Ultrasound Med. Biol, vol. 21, no. 9, pp. 10891100,1995. [71] M. Kinoshita, N. McDannold, F. A. Jolesz, and K. Hynynen, "Noninvasive localized delivery of herceptin to the mouse brain by MRI-guided focused ultrasound-induced blood-brain barrier disruption," in Proc. Natl Acad. Sci., vol. 103, pp. 11719-11723, 2006. [72] F. Y. Yang, W. M. Fu, R. S. Yang, H. C. Liou, K. H. Kang, and W. L. Lin, "Quantitative evalution of focused ultrasound with a contrast agent on blood brain barrier disruption," Ultrasound Med. Biol, vol. 33, pp. 1421-1427, September 2007. 92 [73] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference TimeDomain Method. 3rd ed. Boston: Artech House, 2005. [74] J. Van Bladel, "Electromagnetic fields in the presence of rotating bodies," in Proc. IEEE, vol. 64, pp. 301-318, March 1976. [75] D. E. Lawrence and K. Sarabandi, "Electromagnetic scattering from vibrating penetrable objects using a general class of time-varying sheet boundary conditions," IEEE Trans. Antennas Propag., vol. 54, pp. 2054-2061, July 2006. [76] A. M. Buerkle and K. Sarabandi, "Analysis of acousto-electromagnetic wave interaction using sheet boundary conditions and the finite-difference time-domain method," IEEE Trans. Antennas Propag., vol. 55, pp. 1991-1998, July 2007. [77] A. M. Buerkle and K. Sarabandi, "Analysis of acousto-electromagnetic wave interaction using the finite-difference time-domain method," IEEE Trans. Antennas Propag., vol. 56, pp. 2191-2199, August 2008. [78] Comsol Multiphysics, Comsol Multiphysics Acoustic Module - User Guide. [79] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, Fundamentals of Acoustics. 4th ed. New York: Wiley, 2000. [80] R. F. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGraw-Hill, 1961. [81] A. V. Oppenheim and A. S. Willsky, Signals and Systems. 2nd ed. Englewood Cliffs, NJ: Prentice Hall, 1996. 93 Appendix A: Theoretical Analysis of the Doppler Reconstruction Method A theoretical analysis of the Doppler reconstruction method used in the FDTD-SBC implementation is provided here, with the same 1-D example analyzed as in Section 4.1 and 5.1. Consider that a half space dielectric interface, displacing at an acoustic frequency uim, is illuminated by a continuous plane wave at a microwave frequency UIQ. The scattered microwave field from the interface can be written as Es(z,t) = ft{re-j20ld(t)e+j0lZejuot} (A.l) where F is the reflection coefficient of the dielectric interface, j3\ is the microwave wave number in the medium, and the interface displacement d(t) is d{t) = Amsm(umt) (A.2) The acoustic-modulated scattered field can be expressed as Ea(t)\z=0 = r-cos(w0t - 2/Mm am(umt)) (A.3) If the dielectric interface is unperturbed, then this expression becomes the unperturbed scattered field £b(*)U=o = rcos(u;oi) (A.4) In FDTD-SBC method, the perturbation scattered component is defined as E(t) = Es(t) - E0(t) (A.5) As in the QS FDTD-SBC implementation described in Section 8.2, FDTD-SBC simulations are conducted to compute the perturbation components associated with a series of snapshots at r% (i=l,2,...N) on the acoustic time scale. N is the sampling number within one acoustic period. Then a two-step demodulation process is employed to reconstruct the Doppler component. The computed perturbation component from the FDTD-SBC simulation is a time-domain pulse response which corresponds to a wideband signal in the frequency domain, thus taking the Fourier transform 94 of this perturbation component can extract its phasor component at the center microwave frequency associated with each time snapshot. Next, this phasor component of the perturbation component which varies on the acoustic time scale is demodulated via Fourier transform to obtain the Doppler scattering component. Assume that N snapshots are taken in one acoustic period. The displacement of the dielectric interface associated with each snapshot Tj is 2m d(ri) = Amsm(—),i = l,2,...N (A.6) The total acoustic-modulated scattered field associated with the snapshot Tj is then expressed as 2m Es(t, Ti)\z=0 = rcos(w 0 t - 2f31Amsin(—)), i = 1,2, ...TV (A.7) The perturbation scattered component can then be obtained as ~ E(t,Ti)\z=0 2m = r c o s ( w o * - 2 / M m S i n ( — ) ) - r C os(w 0 t),* = 1,2,...N (A.8) With the assumption that the displacement is much smaller than the microwave wavelength, i.e. Am -C Ai which is true in the proposed application, the expression is then approximated as ~ E{t,Ti)\z=0 2m = Tsm(uj0t)2f31Amsm(—),i = l,2,...N (A.9) The phasor of the perturbation component at the central microwave frequency U>Q associated with each snapshot Tj is extracted via Fourier transform as 27TZ 7T N j E(r i ,w 0 ). = r 2 / 3 i ^ n s i n ( _ ) - - , i = l,2,..JV (A.10) Consider constructing an impulse train to represent the perturbation phasor component that varies on the acoustic time scale N „ L x(t) = r2f31Am(-)J2M^ ) . +oo J2 5^-^T--kT^ (A n - > where Tm = — is the period of the acoustic excitation. This time-varying perturbation phasor component on the acoustic time scale is then demodulated via Fourier transform as X{jku) = T2f31Am(-)J2sin^)uJm J i=l E e fc=-oo ^(^) (AJ2) 95 For the 1 Doppler component, where k = 1 and u = um, the expression above becomes X(jwm) = vr. T—v . ,27ri, _;2*i r2/?!Am(-)o;m^sin(—)e~^ i=l JV/2 iV/2-1 •^ z=l = r2/3^ m (-)o; m V = . sin(-—)e+^] = r2/Mm(-)u;m[J>n(—)e"^+ £ i=0 -2|sin(^)|2 -Tl31Am^)umN J (A.13) Therefore, the magnitude of the demodulated term is \X(jum)\ = r01Am™mN (A.14) Recall that for this canonically deforming case, FDTD-SBC simulation for only one snapshot on the acoustic time scale, i.e. when the dielectric interface is under peak displacement (d(ri) = Am), needs to be conducted to approximate the magnitude of Doppler component. In that case, the computed time-domain perturbation component is E(t,T1)\z=0 = Tsm(uj0t)2p1Am (A.15) Based on which, the magnitude of the 1 s t Doppler component is approximated as |E d (w m )| = r2/3 1 ^ m 7r (A.16) Compare Equations A.14 with A.16. It is evident that the expression obtained as in Equation A.14 calibrated by a normalization factor ^ ^ can be adopted to approximate the magnitude of I s4 Doppler component as in Equation A.16. This factor ^jj normalizes the metric (Equa- tion A.14) by accounting for the chosen total sampling number within one acoustic period in the reconstruction method. 96 Appendix B: List of Authored Publications Peer-Reviewed Publications [1] M. Zhao, S. C. Hagness, B. D. Van Veen, and D. W. van der Weide, "Three-dimensional computational investigation of a focused-acoustic and microwave hybrid modality for breast cancer detection and diagnosis," IEEE Trans. Biomed. Eng., to be submitted. [2] M. Zhao, and S. C. Hagness, "Three-dimensional numerical modeling of electromagnetic scattering from an acoustically excited object using the finite-difference time-domain method," IEEE Trans. Antennas Propag., to be submitted. [3] M. Zhao, J. D. Shea, S. C. Hagness, D. W. van der Weide, and B. D. Van Veen, "Numerical study of microwave scattering in breast tissue via coupled dielectric and elastic contrasts," IEEE Antennas Wireless Propag. Lett., vol. 7, pp. 247-250, 2008. [4] M. Zhao, J. D. Shea, S. C. Hagness, and D. W. van der Weide, "Calibrated free-space microwave measurements with an ultrawideband reflectometer-antenna system," IEEE Microw. Wireless Comp. Lett., vol. 16, no. 12, pp. 675-677, December 2006. [5] M. K. Choi, M. Zhao, S. C. Hagness, and D. W. van der Weide, "Compact mixer-based 1-12 GHz reflectometer," IEEE Microw. Wireless Comp. Lett, vol 15, no. 11, pp. 781-783, November 2005. Conference Presentations [1] M. Zhao*, S. C. Hagness, B. D. Van Veen, and D. W. van der Weide, "Computational study of a focused acoustic and microwave hybrid sensing modality that exploits coupled dielectric and elastic properties contrasts," IEEE APS International Symposium and URSI North American Radio Science Meeting, Charleston, June 1-5, 2009. [2] M. Zhao*, J. D. Shea, S. C. Hagness, B. D. Van Veen, D. W van der Weide, and T. Varghese, "Three-dimensional computational study of microwave scattering in breast tissue via coupled dielectric and elastic contrasts," IEEE AP-S International Symposium and URSI North American Radio Science Meeting, San Diego, July 5 - 1 1 , 2008. 97 [3] S. C. Hagness*, D. W. van der Weide, B. D. Van Veen, M. Zhao, and J. D. Shea"Recent advances in microwave detection of breast cancer: a compact low-cost refiectometer and enhanced sensing modality that couples dielectric and elastic properties contrasts," Department of Defense Breast Cancer Research Program Era of Hope 2008 Meeting, Baltimore, June 25 - 28, 2008. [4] M. Zhao*, J. D. Shea, S. C. Hagness, D. W. van der Weide, B. D. Van Veen, and T. Varghese, "Contrast enhancement for breast cancer detection using acoustic-microwave coupling: a preliminary investigation," URSI North American Radio Science Meeting, Ottawa, Canada, July 22 - 26, 2007 (invited talk). 98 Appendix C: Developed Codes for FDTD-SBC Implementation C.l Determine the surface current locations load iv.mat load jv.mat load kv.mat ib=iv+l; jb=jv+l; kb=kv+l; io=iv/2+l/2; jo=jv/2+l/2; ko=kv/2+l/2; r_phy=5e-3; load dx.mat % dertermine the boundary voxels: numVo=0; iVo= [] ; jVo=[]; kVo= [] ; for ii=iv/2-60:iv/2+60 for jj=jv/2-60:jv/2+60 for kk=kv/2-60:kv/2+60 dr=r_phy-sqrt((ii-io)~2+(jj-jo)~2+(kk-ko)~2)*dx; if dr>0 && dr<dx numVo=numVo+l; iVo(immVo)=ii; jVo(numVo)=jj; kVo(numVo)=kk; end end end end % variables used to determine the surface locations Dexl=r_phy-sqrt((iVo-io).~2+(jVo-0.5-jo).~2+(kVo-0.5-ko). ~2).*dx; Dex2=r_phy-sqrt((iVo-io)."2+(jVo+0.5-jo).~2+(kVo-0.5-ko).~2).*dx; Dex3=r_phy-sqrt((iVo-io).~2+(jVo-0.5-jo).~2+(kVo+0.5-ko)."2).*dx; 99 Dex4=r_phy-sqrt(-(iVo-io) . ~2+(jVo+0.5-jo) .~2+(kVo+0.5-ko)..~2) .*dx; Deyl=r_phy-sqrt Dey2=r_phy-sqrt Dey3=r_phy-sqrt Dey4=r_phy-sqrt iVo-0.5-io iVo+0.5-io iVo-0.5-io iVo+0.5-io ~2+(jVo-jo).~2+(kVo-0.5-ko) ~2+(jVo-jo).~2+(kVo~0.5-ko) ~2+(jVo-jo).~2+(kVo+0.5-ko) ~2+(jVo-jo).~2+(kVo+0.5-ko) "2).*dx; ~2).*dx; "2).*dx; "2).*dx; Dezl=r_phy-sqrt Dez2=r_phy-sqrt Dez3=r_phy-sqrt Dez4=r_phy-sqrt iVo-0.5-io iVo+0.5-io iVo-0.5-io iVo+0.5-io *2+(jVo-0.5-jo).~2+(kVo-ko) *2+(jVo-0.5-jo).~2+(kVo-ko) •2+(jVo+0.5-jo)."2+(kVo-ko) '2+(jVo+0.5-jo).~2+(kVo-ko) "2).*dx; "2).*dx; "2).*dx; ~2).*dx; Dhxl=r_phy-sqrt Dhx2=r_phy-sqrt iVo-0.5-io iVo+0.5-io '2+(jVo-jo) . ~2+(kVo-ko)~2) *dx; '2+(jVo-jo).~2+(kVo-ko) '2) *dx; Dhyl=r_phy-sqrt Dhy2=r_phy-sqrt iVo-io).~2+(jVo-0.5-jo).~2+(kVo-ko) ~2) *dx; iVo-io).~2+(jVo+0.5-jo).~2+(kVo-ko) ~2) *dx; Dhzl=r_phy-sqrt Dhz2=r_phy-sqrt iVo-io).~2+(jVo-jo) '2+(kVo-0.5-ko) *2> *dx; iVo-io).~2+(jVo-jo) *2+(kVo+0.5-ko) '2) *dx; % x-dir interior surface locations (e) iex_intl=iVo(find(Dexl>0)) iex_int2=iVo(find(Dex2>0)) iex_int3=iVo(find(Dex3>0)) iex_int4=iVo(find(Dex4>0)) iex_int=[iex_intl iex_int2 iex_int3 iex_int4]; jex_intl=jVo(find(Dexl>0))-0.5; jex_int2=jVo(find(Dex2>0))+0.5; j ex_int3=jVo(find(Dex3>0))-0.5; j ex_int4=jVo(find(Dex4>0))+0.5; jex_int=[jex_intl jex_int2 jex_int3 jex_int4]; kex_intl=kVo(find(Dexl>0))-0.5; kex_int2=kVo(find(Dex2>0))-0.5; kex_int3=kVo(find(Dex3>0))+0.5; kex_int4=kVo(find(Dex4>0))+0.5; kex_int=[kex_intl kex_int2 kex_int3 kex_int4]; % y_dir interior surface locations (e) iey_intl=iVo(find(Deyl>0))-0.5; iey_int2=iVo(find(Dey2>0))+0.5; iey_int3=iVo(find(Dey3>0))-0.5; iey_int4=iVo(find(Dey4>0))+0.5; iey_int=[iey_intl iey_int2 iey_int3 iey_int4]; jey_intl=jVo(find(Deyl>0)); jey_int2=jVo(find(Dey2>0>); jey_int3=jVo(find(Dey3>0)); jey_int4=jVo(find(Dey4>0)); jey_int=[jey_intl jey_int2 jey_int3 jey_int4]; key_intl=kVo(find(Deyl>0))-0.5; key_int2=kVo(find(Dey2>0))-0.5; key_int3=kVo(find(Dey3>0))+0.5; key_int4=kVo(find(Dey4>0))+0.5; key_int=[key_intl key_int2 key_int3 key_int4]; % z-dir interior surface locations (e) iez_intl=iVo(find(Dezl>0))-0.5; iez_int2=iVo(find(Dez2>0))+0.5; iez_int3=iVo(find(Dez3>0))-0.5; iez_int4=iVo(find(Dez4>0))+0.5; iez_int=[iez_intl iez_int2 iez_int3 iez_int4]; jez_intl=jVo(find(Dezl>0))-0.5; jez_int2=jVo(find(Dez2>0))-0.5; jez_int3=jVo(find(Dez3>0))+0.5; jez_int4=jVo(find(Dez4>0))+0.5; jez_int=[jez_intl jez_int2 jez_int3 jez_int4]; kez_intl=kVo(find(Dezl>0)); kez_int2=kVo(find(Dez2>0)); kez_int3=kVo(find(Dez3>0)); kez_int4=kVo(find(Dez4>0)); kez_int=[kez_intl kez_int2 kez_int3 kez_int4]; '/, x-dir interior surface locations (h) ihx_intl=iVo(find(Dhxl>0))-0.5; ihx_int2=iVo(find(Dhx2>0))+0.5; ihx_int=[ihx_intl ihx_int2]; jhx_intl=jVo(find(Dhxl>0)); jhx_int2=]Vo(find(Dhx2>0)); jhx_int=[jhx_intl jhx_int2]; khx_intl=kVo(find(Dhxl>0)); khx_int2=kVo(find(Dhx2>0)); khx_int=[khx_intl k h x _ i n t 2 ] ; % y - d i r i n t e r i o r surface l o c a t i o n s (h) ihy_intl=iVo(find(Dhyl>0)); ihy_int2=iVo(find(Dhy2>0)); ihy_int=[ihy_intl ihy_int2]; jhy_intl=jVo(find(Dhyl>0))-0.5; jhy_int2=jVo(find(Dhy2>0))+0.5; jhy_int=[jhy_intl jhy_int2]; khy_intl=kVo(find(Dhyl>0)); khy_int2=kVo(find(Dhy2>0)); khy_int=[khy_intl k h y _ i n t 2 ] ; % z - d i r i n t e r i o r surface l o c a t i o n s (h) ihz_intl=iVo(find(Dnzl>0)); ihz_int2=iVo(find(Dhz2>0)); ihz_int=[ihz_intl ihz_int2]; jhz_intl=jVo(find(Dhzl>0)); jhz_int2=jVo(find(Dhz2>0)); jhz_int=[jhz_intl jhz_int2]; khz_intl=kVo(find(Dhzl>0))-0.5; khz_int2=kVo(find(Dhz2>0))+0.5; khz_int=[khz_intl k h z _ i n t 2 ] ; % adjust t o coordinates in FDTD g r i d % x - d i r i n t e r i o r surface l o c a t i o n s (e) 102 iex_int_grid=iex_int; jex_int_grid=jex_int+0.5; kex_int_grid=kex_int+0.5; % y-dir interior surface locations (e) iey_int_grid=iey_int+0.5; jey_int_grid=jey_int; key_int_grid=key_int+0.5; % z-dir interior surface locations (e) iez_int_grid=iez_int+0.5; j ez_int_grid=j ez_int+0.5; kez_int_grid=kez_int; % x-dir interior surface locations (h) ihx_int_grid=ihx_int+0.5; jhx_int_grid=jhx_int; khx_int_grid=khx_int; % y-dir interior surface locations (h) ihy_int_grid=ihy_int; jhy_int_grid=jny_int+0.5; khy_int_grid=khy_int; % z-dir interior surface locations (h) ihz_int_grid=ihz_int; jhz_int_grid=jhz_int; khz_int_grid=khz_int+0.5; % coordinate matrix in FDTD grid °/0 interial surface location (e) cex_int_grid=[iex_int_grid;j ex_int_grid;kex_int_grid] cey_int_grid=[iey_int_grid;jey_int_grid;key_int_grid] cez_int_grid=[iez_int_grid;jez_int_grid;kez_int_grid] % interior surface location (h) chx_int_grid=[ihx_int_grid;jhx_int_grid;khx_int_grid] chy_int_grid=[ihy_int_grid;jhy_int_grid;khy_int_grid] chz_int_grid=[ihz_int_grid;jhz_int_grid;khz_int_grid] C.2 Import the interpolated displacement data from the acoustic-structural simulation °/0 import the Comsol file load fem_app.mat load fa.mat % acoustic frequency wm=2*pi*fa; T=l/fa; M=12; % total snapshots taxis=[0:T/N:T] ; i i = l ; '/, snapshot* ti=taxis(ii); % convert surface location coordinates to % coordinates in acoustic-structural simulations: °/0 x-dir interior surface location (e) xex_int=(iex_int-io).*dx+0.02; yex_int=(jex_int-jo).*dx; zex_int=(kex_int-ko).*dx; % y-dir interior surface location (e) xey_int=(iey_int-io).*dx+0.02; yey_int=(j ey_int-j o).*dx; zey_int=(key_int-ko).*dx; 7, z-dir interior surface location (e) xez_int=(iez_int-io).*dx+0.02; yez_int=(jez_int-jo).*dx; zez_int=(kez_int-ko). *dx; % x-dir interior surface location (h) xhx_int=(ihx_int-io).*dx+0.02; yhx_int=(jhx_int-jo).*dx; zhx_int=(khx_int-ko).*dx; % y-dir interior surface location (h) xhy_int=(ihy_int-io).*dx+0.02; yhy_int=(jhy_int-jo). *dx; zhy_int=(khy_int-ko). *dx; % z-dir interior surface location (h) xhz_int=(ihz_int-io).*dx+0.02; yhz_int=(jhz_int-jo). *dx; zhz_int=(khz_int-ko). *dx; % coordinate matrix for importing data % from acoustic-structural simulations: °/0 for interior surface locations (e) cex_int= [xex_int;yex_int;zex_int]; cey_int=[xey_int;yey_int;zey_int] ; cez_int= [xez_int;yez_int;zez_int]; % for interior surface locations (h) chx_int=[xhx_int;yhx_int;zhx_int]; chy_int=[xhy_int;yhy_int;zhy_int]; chz_int=[xhz_int;yhz_int;zhz_int]; % import interpolated displacement % for x-dir interior surface locations (e) xm_ex_int=postinterp(fem_app,,u_amp_acsld',cex_int); ym_ex_int=postinterp(fem_app,,v_amp_acsld,,cex_int); H- P H- P J J < B o l P 13 13 CD Hi CO M ^_^ ct P H- X O (3* P. P. CO M P- s_^ ^-^ ct s^ ct ^_^ ct P H- p" X 1 «n* *» o *» o CO M O CO M n | pX PJ p" X H- i 1 P H- Hct P P o P- p . h- en1 < P P l 1 1 1 1 13 13 pj Q) QJ B t3- P" 0 B 1 l 1 3 1 3 13 1 1 P P n n P P P P 13 13 B 1 P 13 13 B B 3 1 p 13 13 1 P 13 13 CD CD CD CD Hi Hi /^^ Hi Hi /•—N /-"N /•-N H 13 H. 13 /~N ct ct ct CD CD CD H l-i H 13 1 3 1 3 P H- X 1 CJ* 0 x ct CD P P H- X P 1 1 ct CD P ct P H- P X 1 B N P P p ct ct ct ct II II II II II 1 3 1 3 13 1 3 1 3 O O O O o CO CO CO CO CO c t c t c t ct c t H- M - H - H - H - - P X 1 1 P X i i X 13 <<: 13 p - • ~ N CO p o H- p ct o H- P Ct ct P • - N H- H- << H- H- * CD !->• H- t h-' 00 H- *-•• ct CD ct CD P P ct CD P H- ct CO ct CD P H- ct CO ct CD P H- ct CO Ct CD f-i 13 P H- ct CO H- ^ " N y—\ •X H- ^ - N H- CD o ct P H- \_^ ct P H- ct v p H' P ^^ ct p H- H- ct **s 1 CD i CD p. CO I-1 ct ^_x P H- i CD P. CO M «n "o "n p- CO M s ^-\ CD ^-^ CO p o H- ct M + + -X •—\ /-N H- CD H- X CD ct P H- X 1 H- ct P H- H- ct P H- ct P H- ct P H- ct P H- P ct H- X 1 X CD CD ^ <<: <<: P ct H- CD X 1 X ^-^ ^_1^ X 1 >-L H- P H- P P O CO M P n 1 1P3 1 1 <!»- s~^ Hi CO M - ct s_x P H- X 1 O CD - ^-^ ct p H- ct '—• P H- CD CD X X 1 1 «n «o- ct ^-^ P H- X CD i «n p. M n CO v* E| 1 1 13 P* g 1 13 P 1 n P P v* Hi • ~ \ + P- P. P- 13 13 -X ct P H- 1 CO HM 13 1 3 1 | X 1 ~ E| •X 1 3 ^-N £3- H- ^-N ^ "^ ^ H- 13 •X ct P H- X 1 CD 13 1 + <<S+ N H- H- H P P /~N / " N Hi Hi ^ CD 1 ^- 1 H- CD CD CD CD B B B B 1 1 l l P P P P 1 3 13 1 3 1 3 13 13 13 13 - N CD \ ' H- P ct • •X CD X 13 1 N CD o H- 1 H- ^ - ^ s_^ v_^ *~s ^s \~s *. ' o o H- B ct ct ct ct 1 CD CD CD CD CD H. H. l-i l-i X 13 1 3 13 1 3 + X s; g E| B B B •X -X * ct ct ct ^"N H- /"N -* * •~N H- 13 v_^ 1> > X 13 CD to> X •X P ct H- X 1 1 CD CD + < 1X Ht o> P ct + El 1 CD > B"B to P CD V w / s_^ v—' v_^ **s *^x / ^_x h-' 00 !-»• 00 1 H- ct N^-X P H- 1 n CD h-1 M p . P- - n CD P Hi i-i p CO l-i o ^<! o o H- CD O P- M << << ^ •X- •X • X 1 3 13 13 ct P P ct H- 1 H- <<! 1 1 1 CD CD CD + + X+ N 13 ^ 13 13 H- ^ - N •X- •—s / ^ N H- >-i ct CD P H- CD X P P P P CD c t ri- ri- riCO X X X II ll ll ll 1 1 1 3 1 3 1 3 13 •£> 1 HH - H- O H. O O O ct P P P CO CO CO CO • - v ct ct ct ct ct ct ct N •—N ll P M ^"N P1 I— P ri- CD CD H- ^<! P N N X II 1 3 1 3 1 3 B l-i 1 1 1 CD CD CD CD CD P X X X X H1 1 1 1 I o u . •X- ^ ~ N ) H. H H. H + 1 3 1^ 3 13 1 3 1 3 < *<: *<: <<: i i i / — \ / — \ H i H i H i H i /->. Hi H - H - H - H i /^\ t o> P P P CD CD CD CD CD CD ct ct ct B B B B B B 1 1 1 1 £+ *D C*D 13P 13P 1P3 1P3 13P 13p CD C NJ> X X X 1 3 13 1 3 1 3 1 3 1 3 \_^ 13 13 13 r—\ H - H - H - si < P El < P 1 1 1 1 1 * 13 PJ 1P3 13 P" 0 P 13 1 3 ^ 0^ B^ 1P 1P 1P 13 B 1 1 1 o o P n * c*t CO CO CO n Po Po c*t ct B 1 CD H- ct CO H- ct CO HH P. << 1 < II l-i P. El II M 1 CD X !-»• 00 00 H- P ct P H- 1 N s_^ *^s p ri H- 1 N ^-^ CD H- O l-i Hi 5< 00 N^-X ct P H- I N «o »n* «n• CD CD CD hp. CO p o H- ct P o n h-1 B 1 ct P + + + N ct P ct P to l 1 CD (D p ^ - N H- H- w 3 1 CD CD ^1 l-J- CD N ct M n n CO CO 1 B 13 1 P 13 1 P 1 QJ B P n CD P Hi p CO H O ^ H- CD P HCt H- H- H- CD B N <i <: ^1 *<: •< i 1 •<! 1 ^ 1 ^ N X 1 3 13 1 3 1 1 1 CD CD CD P P P P P P ri- ri- ri- ri- ri- rill ll ll ll ll ll 1 3 1 3 13 13 1 3 1 3 H- O O O O O O y - > N ^-\ X CD ^ CD CD <<! ^1 ^< /^-s II l-i CD P1 h- p 00 (-*• N CD N 1 ct B 1 p 13 13 CD /*-^ Hi CD l-i 13 P ri H- ct i P 13 13 P- p . CO M B l 0 Hi / ^ N 13 H ct CD P H- ct CO f-i ct CO •P ll P ri- H f-i CD P1 h- CD P M <II II E: 00 P H- c t Ct v—/ \~s P H- 13 1 B 1 PJ P 13 13 n P en o h-1 M p . P. CO o B l CD Hi P H- ct CO HH P. N 1 H- ^1 p. 1 CD ^. -^ ~\ I- v_^ ^-x 1 X CD N P ct v—' P H- P. CO M o < SI P 1 1 1 1 3 1 3 13 P* tr P* i P P P B 1 P 13 13 CD Hi P H- ct CO CO H- O M Hi o^ ^^ * ct P H- o o o 00 •X ct P H- 1 3 T) 1 3 H- H- H- * ct P H- + + + H- g B •if £l B •X ct ct / N H- • - N H- P 13 CD B Hi X 1*o 3 13 X) * * H- •X CD X B 1 PJ CD Hi P H- ct CO ct c t ct c t CD CD CD CD H H H H 1 3 13 13 1 3 / — N s~^ ^ " N s~\ P H- ct CO N X 13 13 13 1 1 1 o o n CD CD CD CD CD CD N N N N N N 1 1 1 1 1 H- ct •X B g • - N •X 13 h- 1 O * CD to> X £l N ct P «<! ^ - N ^ X N X N II 1 3 13 B B l-i 1 1 1 l 1 1 CD CD CD CD CD CD CD N N N N P N N 1 1 1 1 p N X P P P P P P CD CD CD c t ri- ri- ri- ri- riN N N II ll ll ll ll ll 1 1 1 1 3 13 13 1 3 1 3 13 H- HH- O O O O O O • ~ \ H CD P> <II + + X+ B B t o> l ^ CD CD CD + < 1N 1N 1N HH- HIO) P p P + c t ct c t P ^ ~ N ct H CO •T! ct II P P> CD O CD P> H Hi P CO H. O I-! H- ct CD P H- HH X 1 p. l-i o ^ II p. E| Hi 1 CD N 1 H- 5-s oca o zp_hx_int=postinterp(fem_app,,w_ph_acsld',chx_int); u=real(xhx_int+xm_hx_int.*exp(i.*(wm*ti+xp_hx_int.*pi./180)))-xhx_int; v=real(yhx_int+ym_hx_int.*exp(i.*(wm*ti+yp_hx_int.*pi./180)))-yhx_int; w=real(zhx_int+zm_hx_iiit.*exp(i.*(wm*ti+zp_hx_int.*pi./180)))-zhx_int; d_hx_int=sqrt(u.~2+v."2+w.~2); °/0 for y-dir interior surface locations (h) xm_hy_int=postinterp(fem_app,'u_amp_acsld',chy_int); ym_hy_int=postinterp(fem_app,'v_amp_acsld',chy_int); zm_hy_int=postinterp(fem_app)'w_amp_acsld',chy_int); xp_hy_int=postinterp(fem_app,'u_ph_acsld',chy_int) yp_hy_int=postinterp(fem_app,'v_ph_acsld',chy_int) zp_hy_int=postinterp(fem_app,'w_ph_acsld',chy_int) u=real(xhy_int+xm_hy_int.*exp(i.*(wm*ti+xp_hy_int.*pi./180)))-xhy_int; v=real(yhy_int+ym_hy_int.*exp(i.*(wm*ti+yp_hy_int.*pi./180)))-yhy_int; w=real(zhy_int+zm_hy_int.*exp(i.*(wm*ti+zp_hy_int.*pi./180)))-zhy_int; d_hy_int=sqrt(u.~2+v.~2+w.~2); % for z-dir interior surface locations (h) xm_hz_int=postinterp(fem_app,,u_a^lp_acsld,,chz_int); ym_,hz_int=postinterp(fem_app, 'v.amp.acsld' ,chz_int) ; zm_hz_int=postinterp(fem_app,,w_amp_acsld',chz_int); xp_hz_int=postinterp(fem_app,'u.ph.acsld',chz_int) yp_hz_int=postinterp(fem_app,'v.ph.acsld',chz_int) zp_hz_int=postinterp(fem_app,'w.ph.acsld',chz_int) u=real(xhz_int+xm_hz_int.*exp(i.*(wm*ti+xp_hz_int.*pi./180)))-xhz_int; v=real(yhz_int+ym_hz_int.*exp(i.*(um*ti+yp_hz_int.*pi./180)))-yhz_int; w=real(zhz_int+zm_hz_int.*exp(i.*(wm*ti+zp_hz_int.*pi./180)))-zhz_int; d_hz_int=sqrt(u.~2+v.~2+w.~2); C.3 Create 'pointers' for checking the neighboring E fields on the same of the boundary % for x-dir interior location (e) n_ex_int=size(iex_int,2); In_ex_int=zeros(n_ex_int,8); for kk=l:n ex int % neighboring Ey if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)-0.5-jo)~2 +(kex_int(kk)-ko)~2)*dx)<r_phy In_ex_int(kk,l)=l; end if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)-0.5-jo)~2 +(kex_int(kk)-ko)"2)*dx)<r_phy In_ex_int(kk,2)=1; end if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)+0.5-jo)~2 +(kex.int(kk)-ko)~2)*dx)<r_phy In_ex_int(kk,3)=1; end if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)+0.5-jo)~2 +(kex_int(kk)-ko)"2)*dx)<r_phy In_ex_int(kk,4)=l; end % neighboring Ez if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)-jo)~2 +(kex_int(kk)+0.5-ko)"2)*dx)<r_phy In_ex_int(kk,5)=1; end if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)-jo)~2 +(kex_int(kk)+0.5-ko)"2)*dx)<r_phy In_ex_int(kk,6)=1; end if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)-jo)~2 +(kex_int(kk)-0.5-ko)~2)*dx)<r_phy In_ex_int(kk,7)=1; end if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)-jo)~2 +(kex_int(kk)-0.5-ko)"2)*dx)<r_phy In_ex_int(kk,8)=1; end end '/, for y-dir interior location (e) ... ... ... ... ... ... ... ... n_ey_int=size(iey_int,2); In_ey_int=zeros(n_ey_int,8); for kk=l:n_ey_int % neighboring Ex if (sqrt((iey_int(kk)-0.5-io)~2+(jey_int(kk)-0.5-jo)~2 +(key_int(kk)-ko)~2)*dx)<r_phy In_ey_int(kk,1)=1; end if (sqrt((iey_int(kk)-0.5-io)~2+(jey_int(kk)+0.5-jo)~2 +(key_int(kk)-ko)~2)*dx)<r_phy In_ey_int(kk,2)=l; end if (sqrt((iey_int(kk)+0.5-io)~2+(jey_int(kk)-0.5-jo)"2 +(key_int(kk)-ko)~2)*dx)<r_phy In_ey_int(kk,3)=l; end if (sqrt((iey_int(kk)+0.5-io)~2+(jey_int(kk)+0.5-jo)~2 +(key_int(kk)-ko)"2)*dx)<r_phy In_ey_int(kk,4)=l; end % neighboring Ez if (sqrt((iey_int(kk)-io)~2+(jey_int(kk)-0.5-jo)~2 +(key_int(kk)+0.5-ko)~2)*dx)<r_phy In_ey_int(kk,5)=1; end if (sqrt((iey_int(kk)-io)"2+(jey_int(kk)+0.5-jo)~2 +(key_int(kk)+0.5-ko)"2)*dx)<r_phy In_ey_int(kk,6)=1; end if (sqrt((iey_int(kk)-io)~2+(jey_int(kk)-0.5-jo)~2 +(key_int(kk)-0.5-ko)"2)*dx)<r_phy In_ey_int(kk,7)=1; end if (sqrt((iey_int(kk)-io)"2+(jey_int(kk)+0.5-jo)"2 +(key_int(kk)-0.5-ko)~2)*dx)<r_phy In_ey_int(kk,8)=l; end end ... ... ... ... ... ... ... ... /o for z-dir interior surface location (e) n_ez_int=size(iez_int,2); In_ez_int=zeros(n_ez_int,8); for kk=l:n_ez_int % neighboring Ex if (sqrt((iez_int(kk)-0.5-io)~2+(jez_int(kk)-jo)~2 +(kez_int(kk)-0.5-ko)~2)*dx)<r_phy In_ez_int(kk,1)=1; end if (sqrt((iez_int(kk)-0.5-io)~2+(jez_int(kk)-jo)~2 +(kez_int(kk)+0.5-ko)~2)*dx)<r_phy In_ez_int(kk,2)=l; end if (sqrt((iez_int(kk)+0.5-io)~2+(jez_int(kk)-jo)~2 +(kez.int(kk)-0.5-ko)~2)*dx)<r_phy In_ez_int(kk,3)=1; end if (sqrt((iez_int(kk)+0.5-io)~2+(jez_int(kk)-jo)~2 +(kez_int(kk)+0.5-ko)"2)*dx)<r_phy In_ez_int(kk,4)=l; end ... ... ... ... % neighboring Ey if (sqrt((iez_int(kk)-io)~2+(jez_int(kk)-0.5-jo)~2 +(kez_int(kk)-0.5-ko)"2)*dx)<r_phy In_ez_int(kk,5)=1; end if (sqrt((iez_int(kk)-io)"2+(jez_int(kk)-0.5-jo)"2 +(kez_int(kk)+0.5-ko)"2)*dx)<r_phy In_ez_int(kk,6)=1; end if (sqrt((iez_int(kk)-io)~2+(jez_int(kk)+0.5-jo)~2 +(kez_int(kk)-0.5-ko)"2)*dx)<r_phy In_ez_int(kk,7)=1; end if (sqrt((iez_int(kk)-io)"2+(jez_int(kk)+0.5-jo)"2 +(kez_int(kk)+0.5-ko)"2)*dx)<r_phy ... ... ... ... In_ez_int(kk,8)=1; end end C.4 Calculate the electric surface currents Js 7„ the electric surface current calculation is done based on the % 3-D FDTD code (Fortran) initially developed by Henri Tandradinata ccc ccc ccc ccc local coordinates of the surface currents: iObsP,jObsP.kObsP displacement imported from the acoustic-structural simulation: delta object radius: a 'pointers' for checking neighbouring E filed: sidDerm do n2=l,obsInPx i=i0bsPx(n2) j=j0bsPx(n2) k=k0bsPx(n2) cnt=sidDermx(n2,l)+sidDermx(n2,2)+sidDermx(n2,3) 1 +sidDermx(n2,4)+sidDermx(n2,5)+sidDermx(n2,6) 2 +sidDermx(n2,7)+sidDermx(n2,8) 1 2 3 4 5 6 7 8 9 imJxl = deltax(n2)*(epsStatM(l)-epsStatM(2)) *(sidDermx(n2Jl)*(Ey(i,j-l,k)-Eyold(i)j-l,k) ) +sidDermx(n2,2)*(Ey(i+l)j-l,k)-Eyold(i+l)j-l,k)) +sidDermx(n2,3)*(Ey(i,j,k)-Eyold(i,j,k)) +sidDermx(n2,4)*(Ey(i+l,j,k)-Eyold(i+l,j,k)) +sidDermx(n2,5)*(Ez(i)j,k)-Ezold(i)j,k)) +sidDermx(n2,6)*(Ez(i+l,j,k)-Ezold(i+l,j,k)) +sidDermx(n2,7)*(Ez(i,j,k-l)-Ezold(i,j,k-l)) +sidDermx(n2,8)*(Ez(i+l,j,k-i)-Ezold(i+l,j,k-l))) /dt/cnt 1 2 imJx2=deltax(n2)*(sigM(l)-sigM(2)) *(sidDermx(n2,l)*Ey(i,j-l,k) +sidDermx(n2,2)*Ey(i+l)j-l)k)+sidDermx(n2,3)*Ey(i,j,k) H- H- H- H- << to C_| B !->• «<! + H- C-i 0 CD M ct ct s_^ h-* h-1 w u H B V_^ (_1. s~\ P / - N *<: B CD i-i a 1 C _ l . • - N v_^ M- 1 "sv w (_1. w s-\ H- c_l. H- / - N * •* m N m N >—• (V <_i. H- /'—v * m X s—' *oo • -J « cn t o w P p to to ^~N << <<: H B v-^ I-"- l_j. /—N •X- v—' ^^ •X- ^^ ^^ * m * m tn m N N X X r~\ H- /-> H- H- M1 * l_i. < !-»• + + I+- ' ** «IV «IV « h-+ w 1 «SV + + \—> '—/ CH-O co H+ cn \n a- O. O Ha P CD CD P. M- ' H- *«IV t_l. u ^-^ \ s to •~\ s W H- h-* | | CO v-^ h-1 • " N (TO 3 H- /•"N * m X V < • < P- P- P- P. O PaD O C CD CD a CD O P H M H H CD B B B B M /•"N • • < B P <<! ^ " N •""N < •< " ^: ^ /—N to /—N P P p P ^ t o tO NO t o P ^-^ * v • *00 * cn <.£. * «t o t o C~ON « CO H- ct H- M- M- M- M- M- H- CO ct P \n to • ™ \ ^< P p to to /—N ^ ^- / p ^_^ s—^ !->• s—' 1 •™N v IV ' s w u ^ - N ^ - N /—N « ^_x V_X v—^ u |_l ^-x i S ^-N |-l v_x V ' v_x l-i ^ ^ IV » C_j. » v_^ ^-^ S ^-\ tO O c+ M CD P, ct p. ^ l«! O M i m cn in x w tv !->• » l_l. - H- CD 1 ct I V ^-^ CD ^-^ i *d - - U . + l-i. 1 I-1 ^-^ HI V t - J . IV H1 " /-N 1 W 1 PI w m N 1 X O ITI o (-• >>! M p . O P. M !->• !V + t—i. (_t. ^-N \ S (_i. s~\ •—N P. (-"• P. H- s_^ 1 < _ l . ^ - N s-^ ^-\ H- ( _ J . <_!. 1 « ** + ^ + (_!. tv1 h-»^-^ (_l. I+ -' I!V IV + «!V ^-^ 1 **IV p. M o s_^ p. (_i. ^ - N ^~N + ^tv h-+ «IV M-1 -IV 1 I - ' ^-^ 1 1 tn 1 N m N o o m N M M I-L (_i. ^~N ^ ^ I-' -X^-^ ^—' * * * * * * ^*> ^* ^^ CD ^~\ ^^ m ^n3 m m m N m m N N N X m x m en ><! ^-s X W ^-\ H- H- H- H- H- /"^ H- H - / ^ c t *~s << P ct CD M P ^~s • - ^ to p P <<: B ^ " N ^ tfs. 00 t o 00 ~J 05 ^cn ^ - N^^- N^^- •™N •—\ to to to s~\ <<! ^ ^<! p P ^ p~> P <<j P. P. P- P. P- P- P. HO p. aD CaD a o CDa CaD CD D C CD CD f-i CD ^ M f-i H ^ ^ B B B B B B B 3 II \p. co + W+ en + w+ w+ C+O ra + ^ *N to II p. * •—\ M- II ! - l B C-l B H- C-l H- «)(XISO)01^UIOK B + OT + C+O Cfl + co M- o H- oo to p. 1=1 p. CD en ^ S~\ H- P. +O C H- p to w * ^ N /•"N P. + CO H- s—> **_• p to ^< B CD a P. CO H- + 00 P to CD p P o t o st. o l-i s* B 00 cn <<! y—N M B ^ <<: M B a CD a CD P. P. t o H- \_/ + C+O w >—' N—S to to S^* T) T) I-* s_>- CO I-' P II t o HII H- O P- v_^ <<: t o 13 -o 1 ^ aCO << *< ^<! h-1 P P /~\ P P CO CO <_1. II II IV C _ i . o O O a- cr cr tv P to ^ - N oD P ^ p—\ C t o t o l-i <-^1 • -rf^ B^ v_/ /—N B l-i ^ <<: B l-i a CD o CD P- P- H- + C+O CO H- <<: B CD t-i o P- CO H- ct II o 1 to I- p P. O P P. CD B X to C-l B + H- M !-»• <-i H- B H- \~^ Ji v_x to t-i B X + + (V IV * N * N P P to to ^i cn a M- + w+ ^^ cn Hp. P. O CD CD •-i Hi B B><! IV IV IV + * m* t*n tn N N ^<! »*_• to P. O CD H- ^-\ P ^-^ P B X M CD 00 0 ) p to •~\ B X CD l-i H- P- P. O O + CO+ CO+ CO H- cn rf=> oo do n2=l,obsInPz i=i0bsPz(n2) j=j0bsPz(n2) k=k0bsPz(n2) cnt=sidDermz(n2,1)+sidDermz(n2,2)+sidDermz(n2,3) 1 +sidDermz(n2,4)+sidDermz(n2,5)+sidDermz(n2,6) 2 +sidDermz(n2,7)+sidDermz(n2,8) 1 2 3 4 5 imJzl = deltaz(n2)*(epsStatM(l)-epsStatM(2)) *(sidDermz(n2,l)*(Ex(i-lJj,k)-Exold(i-l,j,k) ) +sidDermz(n2,2)*(Ex(i-l,j,k+l)-Exold(i-l,j,k+l)) +sidDermz(n2,3)*(Ex(i,j,k)-Exold(i,j,k)) +sidDermz(n2,4)*(Ex(i, j ,k+l)-Exold(i,j,k+l)) +sidDermz(n2,5)*(Ey(i,j-l,k)-Eyold(i,j-l,k)) 6 7 8 9 +sidDermz(n2)6)*(Ey(i)j-l,k+l)-Eyold(i,j-l»^+l)) +sidDermz(n2,7)*(Ey(i)j»k)-Eyold(i)j,k)) +sidDermz(n2)8)*(Ey(i,j)k+l)-Eyold(i,j,k+l))) /dt/cnt 1 2 3 4 imJz2=deltaz(n2)*(sigM(l)-sigM(2)) *(sidDermz(n2,l)*Ex(i-l,j,k) +sidDermz(n2,2)*Ex(i-1,j,k+1)+sidDermz(n2,3)*Ex(i,j,k) +sidDermz(n2,4)*Ex(i)j)k+l)+sidDermz(n2,5)*Ey(i,j-l,k) +sidDermz(n2,6)*Ey(i,j-l,k+l)+sidDermz(n2)7)*Ey(i)j,k) 5 +sidDermz(n2,8)*Ey(i,j,k+1))/cnt imJz=imJzl+imJz2 end do C.5 Determine the index matrix of the neighboring H surface fields % spherical coordinates of the surface field locations % interior surface hx phi_hx_int=acos((ihx_int-io)./sqrt((ihx_int-io)."2 ... +(jhx_int-jo)."2)); % in rad. phi_hx_int(find((ihx_int-jo)<0))=2*pi-phi_hx_int(find((ihx_int-jo)<0)); the_hx_int=acos((khx_int-ko)./sqrt((ihx_int-io). ~2 ... +(jhx_int-jo).~2+(khx_int-ko)."2)); % interior surface hy phi_hy_int=acos((ihy_int-io)./sqrt((ihy_int-io)."2 ... +(jhy_int-jo).~2)); % in rad. phi_hy_int(f ind((jhy_int-j o)<0))=2*pi-phi_hy_int(f ind((j hy_int-jo)<0)); the_hy_int=acos((khy_int-ko)./sqrt((ihy_int-io)."2 ... +(jhy_int-jo).~2+(khy_int-ko)."2)); % interior surface hz phi_hz_int=acos((ihz_int-io)./sqrt((ihz_int-io)."2 ... +(jhz_int-jo).~2)); % in rad. phi_hz_int(find((jhz_int-jo)<0))=2*pi-phi_hz_int(find((jhz_int-jo)<0)); the_hz_int=acos((khz_int-ko)./sqrt((ihz_int-io)."2 ... +(jhz_int-jo).~2+(khz_int-ko).~2)); % determine index matrix of neighboring H surface field: % for interior surface hx n_hx_int=size(ihx_int,2); Ihx_int_nb= [] ; Ihx_int_nb_l=[] ; Ihx_int_nb_2=[] ; for kk=l:n_hx_int de=(ihx_int-ihx_int(kk).*ones(l,n_hx_int))."2+(jhx_int-jhx_int(kk) .*ones(l,n_hx_int)).~2+(khx_int-khx_int(kk).*ones(l,n_hx_int))."2; % avoid choose neighboring locations -> dthe=0 or dphi=0 de(find(the_hx_int==the_hx_int(kk)))=le9; de(find(phi_hx_int==phi_hx_int(kk)))=le9; [Y,I]=sort(de,2); Ihx_int_nb(kk,l)=I(2); for tmp=3:n_hx_int if the_hx_int(I(tmp))~=the_hx_int(I(2)) && ... phi_hx_int(I(tmp))~=phi_hx_int(I(2)) if (the_hx_int(I(tmp))-the_hx_int(kk))*(the_hx_int(I(2)) ... -the_hx_int(kk))<0 && (phi_hx_int(I(tmp))-phi_hx_int(kk)) .. *(phi_hx_int(1(2))-phi_hx_int(kk))<0 Ihx_int_nb(kk)2)=I(tmp); break end end end 13=1(tmp); de=(ihx_int-ihx_int(I(2)).*ones(1,n_hx_int)).~2+(jhx_int-jhx_int(I(2)) .*ones(l,n_hx_int)).~2+(khx_int-khx_int(1(2)).*ones(l,n_hx_int))."2; de(find(the_hx_int==the_hx_int(I(2))))=le9; de(find(phi_hx_int==phi_hx_int(I(2))))=le9; [Y_l,I_l]=sort(de,2); Ihx_int_nb_l(kk,l)=I_l(2); for tmp=3:n_hx_int if the_hx_int(I_l(tmp))~=the_hx_int(I_l(2)) && ... phi_hx_int(I_l(tmp))~=phi_hx_int(I_l(2)) if (the_hx_int(I_l(tmp))-the_hx_int(I(2)))*(the_hx_int(I_l(2)) ... -the_hx_int(I(2)))<0 && (phi_hx_int(I_l(tmp))-phi_hx_int(1(2))) ... *(phi_hx_int(I_l(2))-phi_hx_int(I(2)))<0 Ihx_int_nb_l(kk,2)=I_l(tmp); break end end end de=(ihx_int-ihx_int(I3). *ones(l,n_hx_int)).~2+(jhx_int-jhx_int(I3) .*ones(l,n_hx_int)).~2+(khx_int-khx_int(13).*ones(l,n_hx_int)). ~2; de(find(the_hx_int==the_hx_int(I3)))=le9; de(find(phi_hx_int==phi_hx_int(I3)))=le9; [de_2,I_2]=sort(de,2); Ihx_int_nb_2(kk,l)=I_2(2); for tmp=3:n_hx_int if the_hx_int(I_2(tmp))~=the_hx_int(I_2(2)) && . . . phi_hx_int(I_2(tmp))~=phi_hx_int(I_2(2)) if (the_hx_int(I_2(tmp))-the_hx_int(I3))*(the_hx_int(I_2(2)) ... -the_hx_int(I3))<0 kk (phi_hx_int(I_2(tmp))-phi_hx_int(I3)) ... *(phi_hx_int(I_2(2))-phi_hx_int(I3))<0 Ihx_int_nb_2(kk,2)=I_2(tmp); break end end end end % for interior surface hy n_hy_int=size(ihy_int,2); Ihy_int_nb= []; Ihy_int_nb_l=[]; Ihy_int_nb_2=[]; for kk=l:n_hy_int de=(ihy_int-ihy_int(kk).*ones(l,n_hy_int)).~2+(jhy_int-jhy_int(kk) .*ones(l,n_hy_int)).~2+(khy_int-khy_int(kk).*ones(l,n_hy_int)).~2; % avoid choose neighboring locations -> dthe=0 or dphi=0 de(find(the_hy_int==the_hy_int(kk)))=le9; de(find(phi_hy_int==phi_hy_int(kk)))=le9; [Y,I]=sort(de,2); Ihy_int_nb(kk,l)=I(2); for tmp=3:n_hy_int if the_hy_int(I(tmp))~=the_hy_int(I(2)) kk . . . phi_hy_int(I(tmp))~=phi_hy_int(I(2)) if (the_hy_int(I(tmp))-the_hy_int(kk))*(the_hy_int(I(2)) ... -the_hy_int(kk))<0 kk (phi_hy_int(l(tmp))-phi_hy_int(kk)) ... *(phi_hy_int(I(2))-phi_hy_int(kk)) <0 Ihy_int_nb(kk,2)=I(tmp); break end end end 13=1(tmp); de=(ihy_int-ihy_int(I(2)).*ones(1,n_hy_int)).~2+(jhy_int-jhy_int(I(2)) .*ones(1,n_hy_int))."2+(khy_int-khy_int(I(2)).*ones(1,n_hy_int))."2; de(find(the_hy_int==the_hy_int(I(2))))=le9; de(find(phi_hy_int==phi_hy_int(I(2))))=le9; [Y_l,I_l]=sort(de,2); Ihy_int_nb_l(kk,l)=I_l(2); for tmp=3:n_hy_int if theJhy_int(I_l(tmp))~=the_hy_int(I_l(2)) kk . . . phi_hy_int(I_l(tmp))~=phi_hy_int(I_l(2)) if (the_hy_int(I_l(tmp))-the_hy_int(I(2)))*(the_hy_int(I_l(2)) ... -the_hy_int(I(2)))<0 && (phi_hy_int(I_l(tmp))-phi_hy_int(1(2))) ... *(phi_hy_int(I_l(2))-phi_hy_int(I(2)))<0 Ihy_int_nb_l(kk,2)=I_l(tmp); break end end end de=(ihy_int-ihy_int(I3).*ones(l,n_hy_int)).~2+(jhy_int-jhy_int(I3) ... .*ones(l,ii_hy_int)).~2+(khy_int-khy_int(I3).*ones(l,n_hy_int)).~2; de(find(the_hy_int==the_hy_int(I3)))=le9; de(find(phijiy_int==phi_hy_int(13)))=le9; [de_2,1_2]=sort(de,2); Ihy_int_nb_2(kk,l)=I_2(2); for tmp=3:n_hy_int if the_hy_int(I_2(tmp))~=the_hy_int(I_2(2)) kk . . . phi_hy_int(I_2(tmp))~=phi_hy_int(I_2(2)) if (the_hy_int(I_2(tmp))-the_hy_int(I3))*(the_hy_int(I_2(2)) ... -the_hy_int(I3))<0 && (phi_hy_int(I_2(tmp))-phi_hy_int(I3)) ... *(phi_hy_int(I_2(2))-phi_hy_int(13))<0 Ihy_int_nb_2(kk,2)=I_2(tmp); break end end end end % for interior surface hz ; n_hz_int=size(ihz_int,2); Ihz_int_nb=[] ; Ihz_int_nb_l=[] ; Ihz_int_nb_2=[] ; IpO=[]; Ippi=[]; for kk=l:n_hz_int de=(ihz_int-ihz_int(kk).*ones(l,n_hz_int)).~2+(jhz_int-jhz_int(kk) ... .*ones(l,n_hz_int)).~2+(khz_int-khz_int(kk).*ones(l,n_hz_int))."2; if the_hz_int(kk)==0 I I the_hz_int(kk)==pi ^special location index if the_hz_int(kk)==0 IpO=[IpO kk]; phi_hz_int(kk)=le9; else Ippi=[Ippi kk]; phi_hz_int(kk)=le9; end % avoid choose neighboring locations -> dthe=0 or dphi=0 de(find(the_hz_int==the_hz_int(kk)))=le9; de(find(phi_hz_int==phi_hz_int(kk)))=le9; de(find(the_hz_int==0))=le9; °/„ avoid these two locations de(find(the_hz_int==pi))=le9; [Y,I]=sort(de,2); Ihz_int_nb(kk,l)=I(2); for tmp=3:n_hz_int if abs(phi_hz_int(I(tmp))-phi_hz_int(I(2)))>pi/2 kk abs(phi_hz_int(I(tmp))-phi_hz_int(I(2)))<pi Ihz_int_nb(kk,2)=1(tmp); break end ... end else % avoid choose neighboring locations -> dthe=0 or dphi=0 de(find(the_hz_int==the_hz_int(kk)))=le9; de(find(phi_hz_int==phi_hz_int(kk)))=le9; de(find(the_hz_int==0))=le9; % avoid these two locations de(find(the_hz_int==pi))=le9; [Y,I]=sort(de,2); Ihz_int_nb(kk,l)=K2); for tmp=3:n_hz_int if the_hz_int(I(tmp))~=the_hz_int(I(2)) kk . . . phi_hz_int(I(tmp))~=phi_hz_int(I(2)) if (the_hz_int(I(tmp))-the_hz_int(kk))*(the_hz_int(1(2)) -the_hz_int(kk))<0 kk (phi_hz_int(I(tmp))-phi_hz_int(kk)) *(phi_hz_int(1(2))-phi_hz_int(kk))<0 Ihz_int_nb(kk,2)=I(tmp); break end end end end 13=1(tmp); de=(ihz_int-ihz_int(I(2)).*ones(l,n_hz_int)).~2+(jhz_int-jhz_int(1(2)) .*ones(l,n_hz_int)).~2+(khz_int-khz_int(1(2)).*ones(l,n_hz_int))."2; de(find(the_hz_int==the_hz_int(I(2))))=le9; de(find(phi_hz_int==phi_hz_int(I(2))))=le9; de(find(the_hz_int==0))=le9; % avoid these two locations de(find(the_hz_int==pi))=le9; [Y_l,I_l]=sort(de,2); Ihz_int_nb_l(kk,l)=I_l(2); for tmp=3:n_hz_int if the_hz_int(I_l(tmp))~=the_hz_int(I_l(2)) kk ... phi_hz_int(I_l(tmp))~=phi_hz_int(I_l(2)) if (the_hz_int(I_l(tmp))-the_hz_int(I(2)))*(the_hz_int(I_l(2)) ... -the_hz_int(I(2)))<0 && (phi_hz_int(I_l(tmp))-phi_hz_int(1(2))) .. *(phi_hz_int(I_l(2))-phi_hz_int(I(2))) <0 Ihz_int_nb_l(kk,2)=I_l(tmp); break end end end de=(ihz_int-ihz_int(I3).*ones(l,n_hz_int)).~2+(jhz_int-jhz_int(I3) .*ones(1,n_hz_int))."2+(khz_int-khz_int(13).*ones(1,n_hz_int))."2; de(f ind(the_hz_int==the_hz_int(13)))=le9; de(find(phi_hz_int==phi_hz_int(I3)))=le9; de(find(the_hz_int==0))=le9: % avoid these two locations ... de(find(the_hz_int==pi))=le9; [de_2,1_2] =sort(de,2); Ihz_int_nb=2(kk,l)=I_2(2); for tmp=3:n_hz_int if theJiz_int(I_2(tmp))~=the_hz_int(I_2(2)) ... && phi_hz_int(I_2(tmp))~=phi_hz_int(I_2(2)) if (the_hz_int(I_2(tmp))-the_hz_int(I3))*(the_hz_int(I_2(2)) ... -thejiz_int(13))<0 && (phi_hz_int(I_2(tmp))-phi_hz_int(I3)) ... *(phi_hz_int(I_2(2))-phi_hz_int(13))<0 Ihz_int_nb_2(kk,2)=I_2(tmp); break end end end end % reference: spectial locations: chz_IpO=chz_int_grid(:,IpO); chz_Ippi=chz_int_grid(:,Ippi); Ihz_int_nbpO=Ihz_int_nb(:,IpO); Ihz_int_nbppi=Ihz_int_nb(:,Ippi); Ihz_int_nbpO_l=Ihz_int_nb_l(:,IpO); Ihz_int_nbppi_l=Ihz_int_nb_l(:,Ippi); Ihz_int_nbpO_2=Ihz_int_nb_2(:,IpO); Ihz_int_nbppi_2=Ihz_int_nb_2(:,Ippi); C.6 Calculate the magnetic surface currents Ms % the surface derivative of H field is done based on the % 3-D FDTD code (Fortran) initially developed by Henri Tandradinata ccc ccc ccc ccc local coordinates of the surface currents: iObsP,jObsP,kObsP displacement imported from the acoustic-structural simulation: delta object radius: a index matrix of the neighboring H surface fields: Inb, Inbl, Inb2 ccc for x-component (interior) do n2=l,obsInPx i=i0bsPx(n2) j=j0bsPx(n2) k=k0bsPx(n2) inbla=i0bsPx(Inblx(n2,D) jnbla=j0bsPx(Inblx(n2,l)) knbla=k0bsPx(Inblx(n2,l)) inblb=i0bsPx(Inblx(n2,2)) jnblb=j0bsPx(Inblx(n2,2)) knblb=k0bsPx(Inblx(n2,2)) inb2a=i0bsPx(Inb2x(n2,l)) jnb2a=j0bsPx(Inb2x(n2,l)) knb2a=k0bsPx(Inb2x(n2,1)) inb2b=i0bsPx(Inb2x(n2,2)) jnb2b=j0bsPx(Inb2x(n2,2)) knb2b=k0bsPx(Inb2x(n2,2)) 1 2 d_dsl=l/2*(Hx(inbla,jnbla,knbla)-Hx(inblb,jnblb)knblb)) /(l/2*a*(theMx(Inblx(n2,l))-theMx(Inblx(n2,2))) *(l/2*a*(phiMx(Inblx(n2,l))-phiMx(Inblx(n2,2))))) d_ds2=l/2*(Hx(inb2a,jnb2a,knb2a)-Hx(inb2b)jnb2b,knb2b)) 1 /(l/2*a*(theMx(Inb2x(n2,l))-theMx(Inb2x(n2,2))) 2 *(l/2*a*(phiMx(Inb2x(n2,l))-phiMx(Inb2x(n2,2))))) 1 g_x=l/2*(d_dsl-d_ds2)/(l/2*a*(theMx(Inbx(n2,l)) -theMx(Inbx(n2,2)))*(l/2*a*(phiMx(Inbx(n2,l)) 2 -phiMx(Inbx(n2,2)))))*deltax(n2) end do ccc for y-component (interior) do n2=l,obsInPy i=i0bsPy(n2) j=j0bsPy(n2) 122 k=k0bsPy(n2) inbla=i0bsPy(Inbly(n2,l)) jnbla=j0bsPy(Inbly(n2,l)) knbla=k0bsPy(Inbly(n2,l)) inblb=iObsPy(Inbly(n2,2)) jnblb=j0bsPy(Inbly(n2,2)) knblb=k0bsPy(Inbly(n2,2)) inb2a=i0bsPy(Inb2y(n2,l)) jnb2a=j0bsPy(Inb2y(n2,l)) knb2a=k0bsPy(Inb2y(n2,1)) inb2b=i0bsPy(Inb2y(n2,2)) jnb2b=j0bsPy(Inb2y(n2,2)) knb2b=k0bsPy(Inb2y(n2,2)) 1 2 d_dsl=l/2*(Hy(inblaJnbla,knbla)-Hy(inblb,jnblb,knblb)) /(l/2*a*(theMy(Inbly(n2,l))-theMy(Inbly(n2,2))) *(l/2*a*(phiMy(Inbly(n2,l))-phiMy(Inbly(n2,2))))) 1 2 d_ds2=l/2*(Hy(inb2a)jnb2a)knb2a)-Hy(inb2b)jnb2b)knb2b)) /(l/2*a*(theMy(Inb2y(n2,l))-theMy(Inb2y(n2,2))) *(l/2*a*(phiMy(Inb2y(n2,1))-phiMy(Inb2y(n2,2))))) 1 2 g_y=l/2*(d_dsl-d_ds2)/(l/2*a*(theMy(Inby(n2,l)) -theMy(Inby(n2,2)))*(l/2*a*(phiMy(Inby(n2,1)) -phiMy(Inby(n2,2)))))*deltay(n2) end do ccc for z-component (interior) c expect for special locations: do n2=l,obsInPz i=i0bsPz(n2) j=j0bsPz(n2) k=k0bsPz(n2) inbla=i0bsPz(Inblz(n2,l)> 123 jnbla=j0bsPz(Inblz(n2,l)) knbia=k0bsPz(Inblz(n2,l)) inblb=i0bsPz(Inblz(n2,2)) jnblb=j0bsPz(Inblz(n2,2)) knblb=k0bsPz(Inblz(n2,2)) inb2a=i0bsPz(Inb2z(n2,D) jnb2a=j0bsPz(Inb2z(n2,l)) knb2a=k0bsPz(Inb2z(n2,1)) inb2b=i0bsPz(Inb2z(n2,2)) jnb2b=j0bsPz(Inb2z(n2,2)) knb2b=k0bsPz(Inb2z(n2,2)) 1 2 d_dsl=l/2*(Hz(inbla,jnbla,knbla)-Hz(inblb,jnblb,knblb)) /(l/2*a*(theMz(Inblz(n2,l))-theMz(Inblz(n2,2))) *(l/2*a*(phiMz(Inblz(n2,l))-phiMz(Inblz(n2,2))))) 1 2 d_ds2=l/2*(Hz(inb2a)jnb2a)knb2a)-Hz(inb2b)jnb2b,knb2b)) /(l/2*a*(theMz(Inb2z(n2,l))-theMz(Inb2z(n2,2))) *(l/2*a*(phiMz(Inb2z(n2,l))-phiMz(Inb2z(n2,2))))) 1 2 g_z=l/2*(d_dsl-d_ds2)/(l/2*a*(theMz(Inbz(n2,1)) -theMz(Inbz(n2,2)))*(l/2*a*(phiMz(Inbz(n2,l)) -phiMz(Inbz(n2,2)))))*deltaz(n2) end do c for special locations: do n2=l,obsInPzpO i=i0bsPzp0(n2) j=j0bsPzp0(n2) k=k0bsPzp0(n2) inbla=iObsP(InblzpO(n2,1)) jnbla=j0bsP(Inblzp0(n2,l)) knbla=kObsP(InblzpO(n2,l)) inblb=iObsP(InblzpO(n2,2)) jnblb=jObsP(InblzpO(n2,2)) knblb=k0bsP(Inblzp0(n2,2)) inb2a=i0bsP(Inb2zp0(n2,1)) j nb2a=j ObsP(Inb2zp0(n2,1)) knb2a=k0bsP(Inb2zp0(n2,1)) inb2b=i0bsP(Inb2zp0(n2,2)) jnb2b=jObsP(Inb2zp0(n2,2)) knb2b=k0bsP(Inb2zp0(n2,2)) 1 2 d_dsl=l/2*(Hz(inbla,jnbla,knbla)-Hz(inblb,jnblb,knblb)) /(l/2*a*(theMz(InblzpO(n2,l))-theMz(InblzpO(n2,2))) *(l/2*a*(phiMz(InblzpO(n2,l))-phtMz(InblzpO(n2,2))))) 1 2 d_ds2=l/2*(Hz(inb2a,jnb2a,knb2a)-Hz(inb2b,j nb2b,knb2b)) /(l/2*a*(theMz(Inb2zp0(n2,l))-theMz(Inb2zp0(n2,-2))) *(l/2*a*(phiMz(Inb2zp0(n2,l))-phiMz(Inb2zp0(n2,2))))) g_z=l/2*(d_dsl-d_ds2)/(l/2*(l/2*a*theMz(InbzpO(n2,l)) 1 *(l/2*a*theMz(InbzpO(n2,l))))*abs(sin(phiMz(InbzpO(n2,l)) 2 -phiMz(InbzpO(n2,2)))))*deltaz(n2) end do do n2=l,obsInPzppi i=i0bsPzppi(n2) j=j0bsPzppi(n2) k=kObsPzppi(n2) inbla=iObsP(Inblzppi(n2,1)) jnbla=j ObsP(Inblzppi(n2,1)) knbla=kObsP(Inblzppi(n2,1)) inblb=iObsP(Inblzppi(n2,2)) jnblb=j ObsP(Inblzppi(n2,2)) knblb=kObsP(Inblzppi(n2,2)) inb2a=i0bsP(Inb2zppi(n2,1)) jnb2a=j ObsP(Inb2zppi(n2,1)) knb2a=k0bsP(Inb2zppi(n2,1)) inb2b=i0bsP(Inb2zppi(n2,2)) 125 jnb2b=jObsP(Inb2zppi(n2,2)) knb2b=k0bsP(Inb2zppi(n2,2)) 1 2 d_dsl=l/2*(Hz(inbla,jnbla)knbla)-Hz(inblb)jnblb,knblb)) /(l/2*a*(theMz(Inblzppi(n2,l))-theMz(Inblzppi(n2,2))) *(l/2*a*(phiMz(Inblzppi(n2,1))-phiMz(Inblzppi(n2,2))))) 1 2 d_ds2=l/2*(Hz(inb2aJnb2a,knb2a)-Hz(inb2b,jnb2b,knb2b)) /(l/2*a*(theMz(Inb2zppi(n2,1))-theMz(Inb2zppi(n2,2))) * (l/2*a* (phiMz (Inb2zppi (n2,1)) -phiMz (Inb2zppi (n2,2) ).))) g_z=l/2*(d_dsl-d_ds2)/(1/2*(l/2*a*(pi-theMz(Inbzppi(n2,1))) 1 *(l/2*a*(pi-theMz(Inbzppi(n2,1)))))*abs(sin(phiMz(Inbzppi(n2,1)) 2 -phiMz(Inbzppi(n2,2)))))*deltaz(n2) end do % the magnetic surface current is then calculated in matlab: % parameters load epsl.mat load eps2.mat load sigl.mat load sig2.mat load dx.mat load dt.mat load nL.mat load nR.mat % for x-component (interor) for nn=nL:nR temp=zeros(l,n_hx_int); for mm=l:nn-l f_x=((l/epsi)*exp(-sigl*(dt*(nn-l-mm))/epsl)-(l/eps2) *exp(-sig2*(dt*(nn-l-mm))/eps2)).*ones(1,n_hx_int); Mx(nn,:)=temp-dt.*(mx_l(mm,:)-mx_r(mm,:)).*f_x; temp=Mx(nn,:); end end ... keyboard % for y-component (interor) for nn=nL:nR temp=zeros(l,n_hy_int); for mm=l:nn-l f_y=((l/epsl)*exp(-sigl*(dt*(nn-l-mm))/epsl)-(l/eps2) *exp(-sig2*(dt*(nn-l-mm))/eps2)).*ones(1,n_hy_int); My(nn,:)=temp-dt.*(my_1(mm,:)-my_r(mm,:)).*f_y; temp=My(nn,:); end ... end % for z-component (interor) for nn=nL:nR temp=zeros(l,n_hz_int); for mm=l:nn-l f_z=((l/epsl)*exp(-sigl*(dt*(nn-l-mm))/epsl)-(l/eps2) *exp(-sig2*(dt*(nn-l-mm))/eps2)).*ones(l,n_hz_int); Mz(nn,:)=temp-dt.*(mz_l(mm):)-mz_r(mm,:)).*f_z; temp=Mz(nn,:); end end ...

1/--страниц