Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0001172/00001
## Material Information- Title:
- Monte Carlo Simulation and Experimental Analysis of Dense Medium Light Scattering, with Applications to Corneal Light Scattering
- Creator:
- KIM, KIBUM (
*Author, Primary*) - Copyright Date:
- 2008
## Subjects- Subjects / Keywords:
- Coordinate systems ( jstor )
Cornea ( jstor ) Fiber optic cables ( jstor ) Geometric angles ( jstor ) Lasers ( jstor ) Light scattering ( jstor ) Photons ( jstor ) Simulations ( jstor ) Standard deviation ( jstor ) Wavelengths ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Kibum Kim. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 6/1/2004
- Resource Identifier:
- 53208657 ( OCLC )
## UFDC Membership |

Downloads |

## This item is only available as the following downloads: |

Full Text |

PAGE 1 MONTE CARLO SIMULATION AND EXPERIMENTAL ANALYSIS OF DENSE MEDIUM LIGHT SCATTERING, WITH APPLICATIONS TO CORNEAL LIGHT SCATTERING By KIBUM KIM A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2003 PAGE 2 Copyright 2003 by Kibum Kim PAGE 3 ACKNOWLEDGMENTS I would like to express my deepest gratitude to Dr. David Hahn for his guidance and leadership during this study. He gave me an opportunity to work with him and introduced me to a fascinating area of laser-based diagnostics and bioengineering. Dr. Hahn, a considerable advisor and a mentor, guided me not only about the research but also trivial stuff. I would like to thank my parents for all their unconditional love and support. They trusted me and respected my decisions. I also want to thank my brother for sacrificing everything for me all the time. I would like to express my thanks to Jung, Hyun-Kwon for giving me a piece of advice about debugging and kinematics. Also, I would like to thank Leia Coffey. Although I have never met her before, she left a good reference and suggestions for future work. I am grateful for the concern and advice of Dr. Jorge Carranza, Kenjiro Iida, Brian Fisher, Gregory Yoder, Kathryn Masiello, Vincent Hohreiter, and Allen Ball. Most of all, I would like to thank my wife, Yong-Soon, for her support and encouragement. She gave me a big hand to complete my thesis through proofreading while taking care of all housekeeping, which was deeply appreciated. iii PAGE 4 TABLE OF CONTENTS Page ACKNOWLEDGMENTS.................................................................................................iii LIST OF TABLES.............................................................................................................vi LIST OF FIGURES.........................................................................................................viii ABSTRACT......................................................................................................................xii CHAPTER 1 INTRODUCTION AND BACKGROUND.................................................................1 1.1 The Anatomy of Human Cornea.............................................................................1 1.2 Laser Eye Surgery and Corneal Haze.....................................................................4 1.2.1 Corneal Wound Healing...............................................................................4 1.2.2 Three types of Laser Eye Surgeries..............................................................4 1.2.3 Side Effects of Laser Eye Surgeries.............................................................7 1.2.4 Clinical Assessment of Corneal Haze........................................................10 1.2.5 History of Corneal Haze Measurement Studies.........................................11 1.3 Summary of Light Scattering Theory...................................................................13 1.4 Monte Carlo Simulations......................................................................................19 1.5 Objectives.............................................................................................................20 2 POLYSTYRENE SUSPENSION EXPERIMENTS..................................................22 2.1 The Three-Color Laser System.............................................................................22 2.2 Experimental Method...........................................................................................25 2.2.1 Polystyrene Suspensions............................................................................26 2.2.2 The verification of Suspension Concentration...........................................29 2.2.3 Measurement Method.................................................................................32 2.2.4 The Scattering Sample to PE Reference Ratio...........................................34 2.2.5 Preliminary Testing....................................................................................35 2.3 Data Analysis........................................................................................................35 3 MONTE CARLO SIMULATIONS...........................................................................42 3.1 Development of Overall Algorithm......................................................................42 3.1.1 Probability Distribution of Scattering Angle..............................................43 3.1.2 Probability Distribution of Launch Angle..................................................48 iv PAGE 5 3.1.3 Probability Distribution of Rotational Angle.............................................49 3.1.4 Probability Distribution of Pathlength........................................................51 3.1.5 Absorption..................................................................................................52 3.1.6 Photon Progression.....................................................................................54 3.1.6.1 Coordinate transformation................................................................55 3.1.6.2 The outline of photon progression...................................................58 3.1.7 The Whole Algorithm of the program........................................................59 3.2 Verification of Program........................................................................................61 3.2.1 Verification of the Scattering Angle Distribution......................................61 3.2.2 Verification of the Launch Angle Distribution..........................................62 3.2.3 Verification of the Rotational Angle Distribution......................................63 3.2.4 Verification of the Pathlength Distribution................................................64 3.2.5 Verification of Photon Progression............................................................65 3.3 Monte Carlo Simulation.......................................................................................69 3.3.1 Corneal Simulations...................................................................................69 3.3.2 Polystyrene Suspension Simulations..........................................................78 4 INTEGRATED RESULTS AND DISCUSSION......................................................87 4.1 Integrated Results.................................................................................................87 4.2 Slope of the Two-Color Ratios.............................................................................94 5 CONCLUSIONS AND RECOMENDATIONS.........................................................96 5.1 Conclusions...........................................................................................................96 5.2 Recommendations.................................................................................................99 APPENDIX A CALCULATION OF RADIUS AND CORNEAL EXIT ANGLE..........................101 B FORTRAN MONTE CARLO PROGRAM CODE.................................................102 C PS EXPERIMENTAL DATA CORRECTION........................................................115 LIST OF REFERENCES.................................................................................................119 BIOGRAPHICAL SKETCH...........................................................................................123 v PAGE 6 LIST OF TABLES Table page 1-1 A comparison of the laser eye surgeries.....................................................................6 1-2 Clinical grading scale of corneal haze (0 to 4+) and grading criteria.......................11 1-3 Clinical grading scale of corneal haze (0 to 5+) and grading criteria.......................11 2-1 List of apparatuses for the three-laser system..........................................................26 2-2 Scattering particle density and mean scattering pathlength at 7-days and 3-weeks after laser ablation....................................................................................................26 2-3 Polystyrene microsphere specifications...................................................................28 2-4 Polystyrene suspensions densities for particle size of 7.20 m...............................28 2-5 Mean and standard deviation values of transmission for each suspension..............30 2-6 Values of the dilution factor, the extinction coefficient and the relative standard deviation for each suspension..................................................................................31 2-7 Values of number densities and mean scattering pathlength of suspensions...........32 2-8 The PS signal-to-PE reference ratios for each wavelength and concentrations at each fiber optic. Sigma represents the standard deviation of the signal mean (N=12)......................................................................................................................36 2-9 Two-color ratios for each concentration at each fiber optic. Sigma represents the standard deviation of the signal mean......................................................................39 3-1 The five trials of the mean and standard deviation for sampling the launch angles, along with the average of them................................................................................62 3-2 The five trials of the mean and standard deviation for sampling the pathlength, along with the average of them................................................................................64 3-3 The absolute coordinates of a photon obtained from MC simulation and independent calculation as it progresses through a simulated human cornea..........66 vi PAGE 7 3-4 The relative frequency of 100,000,000 trial photons corresponding to each four case for three wavelength and mean scattering pathlengths.....................................71 3-5 The relative frequency corresponding to collection areas and two color ratios for mean scattering pathlengths.....................................................................................75 3-6 The relative frequency of 10,000,000 photons for four cases at each inner radius..79 3-7 The scattering responses of the PS suspension simulation......................................83 3-8 The two-color ratios of the PS suspension simulation.............................................84 4-1 The correction factor................................................................................................87 4-2 The corrected experimental results (Fiber optic 2)..................................................88 4-3 The corrected experimental results (Fiber optic 3)..................................................88 4-4 The corrected experimental results (Fiber optic 4)..................................................89 C-1 Output of a trial for calculating correction factor...................................................117 vii PAGE 8 LIST OF FIGURES Figure page 1-1 Schematic of the eye..................................................................................................1 1-2 Schematic of microscopic composition of the human cornea....................................3 1-3 Light scattering response to an incident light...........................................................14 1-4 Mie light scattering particle geometry.....................................................................16 2-1 Picture of experimental setup of the three-color laser system.................................23 2-2 Schematic of experimental setup of the three-color laser system............................24 2-3 Diagram of the fiber optic........................................................................................24 2-4 The placement of fiber optics in the aluminum bar..................................................25 2-5 Expression of scattering particles in suspension as a volume fraction of one.........27 2-6 Schematic of transmission measurement.................................................................29 2-7 Pictures of measurements for PE reference (Left) and PS suspension (right).........33 2-8 A representative scattering spectrum of PE measurement.......................................33 2-9 The absolute peak intensity......................................................................................35 2-10 The scattering signal-to-reference ratios at fiber optic 2. The error bars denote standard deviations12 sets of experimental mean performed on different days....37 2-11 The scattering signal-to-reference ratios at fiber optic 3. The error bars denote standard deviations12 sets of experimental mean performed on different days....37 2-12 The scattering signal-to-reference ratios at fiber optic 4. The error bars denote standard deviations12 sets of experimental mean performed on different days....38 2-13 The two color ratios at fiber optic 2. The error bars denote standard deviations12 sets of experimental mean performed on different days..........................................39 viii PAGE 9 2-14 The two color ratios at fiber optic 3. The error bars denote standard deviations12 sets of experimental mean performed on different days..........................................40 2-15 The two color ratios at fiber optic 4. The error bars denote standard deviations12 sets of experimental mean performed on different days..........................................40 3-1 Diagram of a light scattering from a scattering particle...........................................43 3-2 Differential scattering cross section for three wavelengths......................................45 3-3 Flow-Chart how to set up probability array for scattering angle theta. ...................46 3-4 Distribution of the laser light exiting from the fiber optic probe tip........................48 3-5 Representation of the rotational angle......................................................................50 3-6 Transmission of the human cornea...........................................................................53 3-7 Simple geometry of the simulated human cornea....................................................55 3-8 Depiction of point transformation............................................................................56 3-9 Flow chart of Monte Carlo program algorithm (Sample array generation).............59 3-10 Flow chart of Monte Carlo program algorithm (Photon progression through the cornea)......................................................................................................................60 3-11 Verification of the scattering angle theta.................................................................61 3-12 Verification of sampling algorithm for launch angle alpha. Two sets of MC trials of 65,000 samples with the Gaussian probability function..........................................63 3-13 Verification of sampling algorithm for rotational angle beta. Two sets of MC trials of 65,000 samples with the uniform probability function........................................64 3-14 Verification of sampling algorithm for pathlength. Two sets of MC trials of 65,000 samples each with the lognormal probability function............................................65 3-15 Tracks of four back-exiting photons........................................................................67 3-16 Tracks of four front-exiting photons........................................................................68 3-17 Simulation of photons absorbed...............................................................................68 3-18 Two sets of program output for 670 nm wavelength and 15.0 m mean scattering pathlength based on 100,000,000 trials each............................................................70 ix PAGE 10 3-19 Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 5.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength...............................................................................................................72 3-20 Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 10.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength...............................................................................................................72 3-21 Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 15.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength...............................................................................................................73 3-22 Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 20.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength...............................................................................................................73 3-23 Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 25.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength...............................................................................................................74 3-24 Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 30.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength...............................................................................................................74 3-25 The two color ratio for each mean scattering pathlength at fiber optic 2.................76 3-26 The two color ratio for each mean scattering pathlength at fiber optic 3.................76 3-27 The two color ratio for each mean scattering pathlength at fiber optic 4.................77 3-28 The relative frequency of photon exiting through front surface of 10,000,000 trials at each inner radius...................................................................................................79 3-29 The front-exiting angle distribution of C4 suspension for three wavelengths. Mean scattering pathlength corresponding to C4 suspension is 120.3 m........................80 3-30 The front-exiting angle distribution of C3 suspension for three wavelengths. Mean scattering pathlength corresponding to C3 suspension is 59 m.............................81 3-31 The front-exiting angle distribution of C2 suspension for three wavelengths. Mean scattering pathlength corresponding to C2 suspension is 34.4 m..........................81 3-32 The front-exiting angle distribution of C1 suspension for three wavelengths. Mean scattering pathlength corresponding to C1 suspension is 27.2 m..........................82 3-33 The two color ratios at fiber optic 2.........................................................................84 3-34 The two color ratios at fiber optic 3.........................................................................85 x PAGE 11 3-35 The two color ratios at fiber optic 4.........................................................................85 4-1 The two-color ratio values (532 nm/594 nm) at fiber optic 2 as a function of concentration for the MC simulation and the PS suspension experiment................89 4-2 The two-color ratio values (594 nm/670 nm) at fiber optic 2 as a function of concentration for the MC simulation and the PS suspension experiment................90 4-3 The two-color ratio values (532 nm/670 nm) at fiber optic 2 as a function of concentration for the MC simulation and the PS suspension experiment................90 4-4 The two-color ratio values (532 nm/594 nm) at fiber optic 3 as a function of concentration for the MC simulation and the PS suspension experiment................91 4-5 The two-color ratio values (594 nm/670 nm) at fiber optic 3 as a function of concentration for the MC simulation and the PS suspension experiment................91 4-6 The two-color ratio values (532 nm/670 nm) at fiber optic 3 as a function of concentration for the MC simulation and the PS suspension experiment................92 4-7 The two-color ratio values (532 nm/594 nm) at fiber optic 4 as a function of concentration for the MC simulation and the PS suspension experiment................92 4-8 The two-color ratio values (594 nm/670 nm) at fiber optic 4 as a function of concentration for the MC simulation and the PS suspension experiment................93 4-9 The two-color ratio values (532 nm/670 nm) at fiber optic 4 as a function of concentration for the MC simulation and the PS suspension experiment................93 4-10 The slope values of the two-color ratios as function of fiber optic position corresponding to collection angle in the PS suspension experiment.......................95 4-11 The slope values of the two-color ratios as function of fiber optic position corresponding to collection angle in the PS suspension simulation.........................95 A-1 Geometry of radius and angle ..............................................................................101 xi PAGE 12 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science MONTE CARLO SIMULATION AND EXPERIMENTAL ANALYSIS OF DENSE MEDIUM LIGHT SCATTERING, WITH APPLICATIONS TO CORNEAL LIGHT SCATTERING By Kibum Kim August 2003 Chair: David Hahn Major Department: Mechanical and Aerospace Engineering Recently, the refractive laser surgery has been very common for those who want to correct visual refraction. However, the surgery potentially has a negative effect which is the presence haze of the cornea. Such corneal haze, while often transient, may lead to the loss of corneal transparency after surgery. Sometimes, corneal haze could also result from the diseases associated with the cornea. Despite its clinical significance, there is no effective method to quantify haze objectively. For the benefit of those who suffer from corneal haze, it is necessary to develop a technique that could unveil ambiguous points about corneal haze. Such a technique could be more valuable to understand the processes of tissue wound healing and drug testing as well as to enable ophthalmologists to better monitor the progression of haze. Only a few techniques of clinical haze assessment such as high frequency ultrasound, Scheimpflug photography and light scattering measurements have so far been reported; however, none of these methods are widely applied to the clinic, and corneal xii PAGE 13 haze assessment post-photorefractive keratectomy (PRK) still is dependent primarily on subjective examination. To more accurately recognize the pathogenic mechanisms of corneal haze, there is a basic need for an objective in vivo method to quantitatively evaluate corneal haze. As a objective technique to determine corneal haze with an absolute measurement, this study proposes to use a multi-color light scattering signal ratio for the purpose of minimizing individual-dependent scattering variations and for providing an accurate, precise, quantitative measurement of corneal haze that does not require a pre-operative or average normal corneal reference. This study was conducted in terms of the current view that corneal light scattering (haze) is produced by fibroblastic keratocytes. Experiments with polystyrene (PS) suspensions, which simulated the dense medium scattering environment of the human eye, were also performed and analyzed. A two-color scattering ratio, a ratio of the scattering responses of different wavelengths, was calculated and it was found that the ratio increases as scattering particle density increases. Finally, Monte Carlo simulations were performed to verify the PS suspension results and to comprehend the complex multiple-scattering problem, thereby providing the framework for future optical probe design, including optimization studies. Overall, the multi-color scattering system is a promising tool to quantify corneal haze. The results indicate that the two-color normalization of the scattering signal is successful in producing a reasonably precise metric of corneal haze and minimizing person-to-person variations. xiii PAGE 14 CHAPTER 1 INTRODUCTION AND BACKGROUND 1.1 The Anatomy of Human Cornea The cornea, a dome-shaped structure covering both the pupil and the iris, supports the eye mechanically and protects it chemically from outer stimulation. The external tunic of the eye consists of one-sixth of the transparent cornea at the front and five-sixths of the opaque sclera at the limbus. A schematic of the eye is shown in Figure 1-1. Sclera Iris Cornea Lens Figure 1-1. Schematic of the eye. Unlike most other tissues in the body, the cornea does not contain any blood vessels because it has a refractive function. In other words, it has become an avascular structure to perform the optical function efficiently. Therefore, the cornea assimilates required oxygen from the atmosphere and absorbs nutrition from the aqueous humor directly, while maintaining its front surface in a clean and smooth condition by the 1 PAGE 15 2 blinking action and tears from lacrimal glands. The cornea is the major refractive surface along with the lens in the eye, and the front corneal surface is endowed with a refractive index of approximately 1.337 and a power of 48 dioptiers (Klyce & Beuerman, 1998; Martini, 2001; Niemz, 1996). The rear surface of the cornea is almost spherical whereas the front surface of it is not uniformly spherical as shown in Figure 1-1. The shape of the anterior corneal surface is oval; thus, not only is the cornea asymmetric but also the thickness of it is not constant. The average thickness is 520 m centrally and 650 m peripherally or more. The average horizontal diameter is 12.6 mm, and the average vertical diameter is 11.7 mm. The outer radius of curvature is approximately 8 mm, as estimated at the center of the front surface (Klyce & Beuerman, 1998; Niemz, 1996). Microscopically, the corneal tissue is composed of five basic layers: the epithelium, Bowmanâ€™s layer, stroma, Descemetâ€™s membrane, and the endothelium. The epithelium, constituting approximately 10 percent of the corneal thickness, is the most anterior membrane which functions as a barrier between the cornea and the environment. The cornea is extremely sensitive to external shock because numerous microscopic nerve endings compose the epithelium. Bowmanâ€™s layer is made up of strongly packed protein collagen fibers. It is situated along the surface of the cornea and contributes morphologically to the corneal integrity. Once Bowmanâ€™s layer is injured, it cannot be restored and can be scarred, which is possibly a cause of loss the vision ability. Under Bowmanâ€™s layer is the stroma, comprising approximately 90 percent of the corneal thickness, which consists of 78 percent of water and 16 percent of lamellar collagen fibers. The collagen matrix has a unique size and arrangement that provides enough PAGE 16 3 mechanical strength and transparency for the cornea. Descemetâ€™s membrane, a layer between the stroma and the endothelium, is composed of a thin but strong sheet of collagen fiber tissues. This is created recurrently by the endothelial cells and able to be regenerated. The endothelium, the posterior layer, plays an important role in maintaining corneal dehydration. Because the surplus fluid into the stroma can become a cause of the corneal haze and loss of transparency, the endothelium conducts the pumping mechanism so that the cornea can keep a balance of fluid volume and maintain its transparency. Endothelial cells are not able to be regenerated if damaged. The microscopic composition of the human cornea is shown in Figure 1-2. Bowmanâ€™s layer (15m) Stroma (450m) Endothelium Epithelium (40~50 m) Descemetâ€™s membrane Figure 1-2. Schematic of microscopic composition of the human cornea. PAGE 17 4 1.2 Laser Eye Surgery and Corneal Haze 1.2.1 Corneal Wound Healing The cornea is likely to be affected by impairment that results from physical or chemical impacts under any circumstances. Additionally, epithelial defects, epithelial and superficial stromal defects, deep stromal defects and full-thickness defects are mainly the four types of the common corneal trauma. In particular, deep stromal defects can follow after laser refractive surgery, which is one of the potential side effects. As for the corneal healing process, it is influenced by not only the types of trauma mentioned above, but also diverse elements like age, nutrition, trauma cause, infection, and vascularization. Most of all, the keratocyte reaction to trauma plays a vital role in the efficiency of the corneal healing process. In the process of keratocyte response to trauma, stromal keratocytes within the impaired portion of the stroma are activated and turned into fibroblasts, which changes their shape. Consequently, the nucleoli of stromal keratocytes get developed larger than their original size. These activated keratocytes yield collagenous scar tissue during the duration of deep stromal healing. Keratocyte reaction to trauma consequently affects scar formation and the final outcome of healing. Furthermore, keratocytes contribute to yielding fibronectin that is necessary for rebuilding up the corneal epithelial (Parrish & Chandler, 1998). 1.2.2 Three types of Laser Eye Surgeries Methods of refractive vision correction have been developed and have become very popular since keratomileusis was originated in the 1950s. At the initial times, a tiny diamond blade was used to perform the ablation of tissue. Photo-refractive keratectomy (PRK) and Laser in situ keratomileusis (LASIK) were practiced after an excimer laser was introduced. Recently, Laser epithelial keratomileusis (LASEK) has moved into the PAGE 18 5 spotlight as a more advanced technique. These techniques such as PRK, LASIK and LASEK use the excimer laser to reshape the cornea about the optical axis. The revised cornea obtains the correct geometry, which leads to corrections in the refractive properties of the cornea. The excimer laser used in refractive surgeries is a 193 nm ultraviolet, pulsed (5 to 60 Hz) argon-fluoride laser. It has the photon energy of approximately 6.4 eV, which is greater than the bond energy of the collagen protein matrix of the eye. In the area exposed to the laser, the bonds are easily photo-degraded by the laser pulse and the tissues are detached without thermal damage to the remaining tissue. This is known as photodecomposition, and is sometimes termed cold ablation because there is no thermal effect in using the excimer laser (McDonald & Chitkara, 1998; Manche et al., 1998). In a study led by Eiferman (1991), it was shown that the excimer laser does not result in any significant side effects on the eye. Also, the corneal absorption coefficient for 193 nm radiation is so high (~40,000 cm-1) that the incident pulse is absorbed within a cellular layer of the ablated tissue. Although the basic principles of laser eye surgeries are the same, their applications are different in the process of medical practice. For the first time, PRK was introduced in the beginning of 1980s and was devised by Trokel. In PRK applications, the epithelium over the ablation area is eliminated chemically (20% alcohol), mechanically (scraped or brushed away) or by the laser itself, and then the exposed corneal tissue is ablated by the laser. The removed epithelium regenerates naturally. In 1990, LASIK was devised by Pallikaris as an advanced technique that combines keratomileusis with PRK. In the LASIK application, the cornea is anesthetized and cut by PAGE 19 6 a microkeratome to make an anterior corneal flap (thickness of 160~180 m) that includes the epithelium, Bowmanâ€™s layer and a part of stromal tissue. The flap is moved aside, and the stromal corneal tissue is ablated by the laser. Finally, the corneal flap is replaced (Manche et al., 1998). LASEK was developed by Camellin in 1990 as a novel method to take some of the advantages of both PRK and LASIK. The application is the same as LASIK except that an anterior corneal flap (thickness of 60~70 m), which only includes the epithelium, is made by Epi-trephine or using a dilute alcohol solution, which has brought a decrease in serious complications and risks involved with the corneal flap in LASIK. Also, it is different from the PRK in that the epithelium is preserved and replaced after the procedure, which results in a more painless recovery and a faster return of visual acuity. Considering that the occurrence of corneal haze is somewhat high compared to LASIK according to some studies, LASEK is more appropriate in the case that the cornea is too thin to apply LASIK. Table 1-1 shows a comparison of the laser eye surgeries. Table 1-1. A comparison of the laser eye surgeries. PRK LASIK LASEK Objective (dioptier) 3.0 ~ -4.0 -6.0 ~ -15.0 -3.0 Pain Severe Slight Moderate Procedure Simple Complicated Complicated Surgery time < 10Min. ~ 15Min. ~ 15Min. Healing time About 3 months 1 ~ 4 weeks 1 month Stabilization of Refraction 1 ~ 4 weeks 1 Week 1 ~ 2 Weeks Experience ~ 15years Recent Recent FDA recognition Yes (1995) Not yet Not yet Cost Inexpensive Expensive Expensive Damage of Bowmanâ€™s layer Much Little Much Possibility of Side Effect Myopia regression Much Little Little Corneal haze Much Little Much Keratectasia Little Much Little PAGE 20 7 1.2.3 Side Effects of Laser Eye Surgeries There are several complications of laser refractive procedures: refractive regression (myopia regression), corneal haze, glare, halo and decreased contrast following laser eye surgery. Refractive regression occurs when the revised cornea re-thickens gradually, thereby losing 0.5 to 3.0 dioptiers (D) from the corrected refractive index within the first year after the surgery. A higher initial correction of corneal refraction is inclined to bring more myopic regression. It often takes a long time to stabilize post-operative refraction of the cornea due to myopic regression slowly occurring over the period of 3 or 5 years (Moller-Pedersen et al., 2000). Typically, myopic regression is more likely to be developed in the wave of PRK application due to much more damage of epithelial tissue, compared to other therapies. As a result of myopic regression, it is hard to attain stabilization of refraction completely through the first PRK; therefore, a follow-up surgery is inevitable. Another frequently appearing side effect is corneal haze, which is the optic disturbance in corneal transparency consequently following surgery. This negative blur of the cornea appears within the first four months post-operative and usually dissipates within the first 12 months post-operative. However, corneal haze may persist in roughly 5% of all PRK patients (Moller-Pedersen et al., 2000). Many researchers have studied to figure out the cause of corneal haze related to refractive surgery since the incidence of the haze was first notified in 1988. Seiler et al. reported that human corneal haze had trivial grading levels of haze and disappeared by the sixth month following surgery (Seiler et al., 1990). Another study of 120 subjects exhibited a different trend, namely that 92% showed some grade of stromal haze beginning at 4~6 weeks after operation and peaking at about 6 months after it (Gartry et PAGE 21 8 al., 1991). Similarly, corneal haze starting six months post-operative was observed in 63% of myopia-corrected eyes where 54% of the cases were trace, 7.3% were mild, and 2.3% were moderate. In the similar study, 43.9% had trace haze one year post-operation (Odrich & Greenberg, 1998). In a study conducted using rabbit eyes, haze started two weeks post-operative and reached a maximum at 4-weeks post-operative (Kaji et al., 1998). Theories that back scattering is the cause of corneal haze have been generally accepted and changed over the years. Keratocytes in the middle of wound healing become fibroblastic and then begin to yield collagen a couple of days post-operative (McDonald & Chitkara, 1998). Furthermore, a new theory supports that corneal haze is related to enhancement of the back scattered light in the ablated area of keratocytes rather than changes in stromal extracellular fibrosis (Moller-Pedersen et al., 2000). Currently, the activated keratocytes (fibroblastic keratocytes or myofibroblasts) are regarded to play a major role as light scattering centers. For the purpose of assuring this theory of keratocyte scatters, Moller-Pedersen et al. and Jester et al. used confocal microscopy through focusing (CMTF) for myopic patients. Under the normal pre-surgery condition, images of the corneal stroma showed the transparent cell bodies of keratocytes and no scattered light. On the other hand, those images showed reflected light from the nuclei of corneal keratocytes after PRK. In addition, greater haze regions corresponded with regions where more back scattering occurs in keratocytes (Jester et al., 1999a). A later study conducted using CMTF for 17 myopic patients proved that higher numbers of wound healing keratocytes resulted in an increase in cellular reflectivity. The study also confirmed that corneal haze originated from activated fibroblasts, and it provided the PAGE 22 9 evidence that the densities of fibroblastic keratocytes affected the temporal appearance and disappearance of the haze (Moller-Pedersen et al., 2000). Del Val et al. reported that the diameter of these light-scattering keratocytes ranged between 4 and 10 m and had a relative refractive indices ranging from1.01 to 1.08 (del Val et al., 2001). In other studies, the information on the size was verified by recorded images (Jester et al., 1999a, Jester et al., 1999b). Linear scanning densitometry was used to study the light scattered through the cornea by Smith et al. They observed that the light scattering was not uniform, but distributed dissipatedly over the cornea; in other words, it was concentrated at the anterior and posterior rather than at the center of the cornea. Moreover, light scattering was distributed more intensely at the anterior than at the posterior stroma because collagen fibres are more irregularly regenerated at the anterior than at the posterior stroma. It was also noted that most scattering occurs at the interface of two mediums if their refractive indices have a large mismatch. It could be taken as an example that light scattering is concentrated at the epithelial surfaces. On the other hand, some ophthalmologists have treated corneal haze with corticosteroids and reported that corneal haze resulting from back scattered light might not be linked directly to the loss of visual sensitivity. For instance, 93% of patients did not demonstrate that corneal haze caused a loss of their sensitivity (Eiferman et al., 1991; Gartry et al., 1991; Moller-Pedersen et al., 2000). Until now, structural changes in the stroma of the cornea have been emphasized as a first-hand cause of corneal haze in many studies. However, the damage of Bowmanâ€™s layer may also be a factor of corneal haze, which can be corroborated by the comparison PAGE 23 10 of the refractive surgeries shown in Table 1-1 (in section 1.2.2). Table 1-1 indicates that LASIK has a low incidence of corneal haze compared to others because Bowmanâ€™s layer is rarely affected in LASIK. In general, Bowmanâ€™s layer remains scarred and is unable to be regenerated after ablation, while the epithelium is able to be quickly restored and almost completed within the first week post-operation. Therefore, it is suspected that the damage of Bowmanâ€™s layer may result in the corneal haze regardless of the finding that Bowmanâ€™s layer is not necessary for re-epithelialization, and the subepithelial haze appears to have little effect on visual acuity or contrast sensitivity according to a study of six myopic patients conducted by Eiferman and co-workers. From the other prospects, corneal haze can also result from corneal diseases such as herpetic keratitis, corneal ulceration, corneal dystrophies, and other types of corneal trauma. This corneal haze can obscure patientsâ€™ vision, and even cause a loss of visual ability. Therefore, such haze is important even though only a small percentage of patients suffer from such diseases. In short, it is significant to do further studies regarding the factors and healing processes of corneal haze. Such studies can be aided by the ability to quantify corneal haze. 1.2.4 Clinical Assessment of Corneal Haze The standard methods to evaluate corneal haze have been developed by ophthalmologists. Some ophthalmologists employ a slit-lamp microscope to determine and grade the haze of the cornea on the scale from 0 to 4+, as shown in Table 1-2 (del Val et al., 2001; Fantes et al., 1990). Meanwhile, others use a scale from 0 to 5+ as listed in Table 1-3 (Andrade et al., 1990). It is also prevalent to use Food and Drug Administration standardized categories: PAGE 24 11 grade 0, clear; grade 1, minimal haze; grade 2, trace haze; grade 3, mild haze; grade 4, moderate haze; grade 5, marked haze). Table 1-2. Clinical grading scale of corneal haze (0 to 4+) and grading criteria. Clinical Scale Grading Criteria 0 The cornea is completely clear 0.5+ Faint, trace haze is visible by oblique illumination 1+ Minimal haze is visible with direct diffuse illumination 2+ Mild haze is visible with direct focal-slit illumination 3+ Moderate haze is visible which partially obscures the iris details 4+ Severe haze is visible which completely obscures the iris details Table 1-3. Clinical grading scale of corneal haze (0 to 5+) and grading criteria. Clinical Scale Grading Criteria 0 Clear, possibly with faint haze 1+ Mild reticulated haze not affecting refraction 2+ Moderate haze, refraction possible but difficult 3+ Opacity prevents refraction; anterior chamber easily viewed 4+ Opacity impairs view of anterior chamber 5+ Unable to see anterior chamber Not only are these kinds of methods subjective, but repeatability is also unreliable. As mentioned, clinical assessment of corneal haze has not been objective in the absence of the reliable method; therefore, this study is primarily aimed to design an objective measurement technique to quantify corneal haze. 1.2.5 History of Corneal Haze Measurement Studies Andrade et al. attempted to measure the opacity of the human cornea after PRK by using an opacity lensometer that was used for measuring the opacity of the lens. In the study, it was supposed that corneal haze would be a cause of the opacity changes during the post-operative period if the opacity values of the lens remained constant. However, PAGE 25 12 this study was not much helpful for correlating the results to clinical grades due to large variations observed from eye to eye (Andrade et al., 1990). Lohmann et al. utilized a slit-lamp light source and a CCD camera for quantifying human corneal haze. Reflected light was distinguished from scattered light by using polarizers. This study also could not reveal objective values because of variations from individual to individual and large standard deviations caused by a failure to completely account for surface reflection (Lohmann et al., 1992). A scatterometer consisting of a modified slit-lamp microscope, a fiber optic, filters for wavelength selection, and a photomultiplier detector has also been employed as a measuring equipment for corneal light scattering. Jain et al. used it to measure back-scattered light from pigmented rabbit eyes and to analyze how the corneal haze responses to LASIK and PRK. An average of scattering obtained from 100 normal rabbit eyes was used to normalize the degree of light scattering measured after operation. This studyâ€™s finding supported the hypothesis that PRK induces much more corneal scatter than LASIK. Nevertheless, the method of this study was inadequate to be used as a clinical tool because it did not suggest any adoptable clinical grade, which was attributable to the large standard deviations in the scattering index (Jain et al., 1995). With the scatterometer, Braunstein et al. studied that a scattering index comparing the scattering of an ablated cornea to an average scattering of normal corneas could be used as an appropriate standard for objective scattering measurement. However, the scattering index turned out not to correlate well with clinical haze grading (Braunstein et al., 1996). PAGE 26 13 Using confocal microscopy through focusing (CMTF) contributed to studying the human corneal haze despite the shortcomings of expense and precision. CMTF enabled the thickness of the epithelium and stroma to be analyzed, and also made a difference in classifying clinical haze values as well as levels of subclinical haze. Besides, it provided further scientific evidence that deeper ablation led to higher grades of corneal haze. It also proved that the activation of keratocytes occurred more at the anterior 100 m than at the posterior in the cornea (Moller-Pedersen et al., 2000). Finally, it demonstrated that light scatters at the air-cornea and the water-cornea interface. For the purpose of measuring light scattering, Del Val et al. found a relation between transmittance and keratocyte density by applying transmission measurements to Iber Braun hens. In using these transmission measurements, the cornea was to be removed from the original position; therefore, this method does not fit in haze measurements for the human cornea under clinical conditions. (del Val et al., 2001). 1.3 Summary of Light Scattering Theory Light scattering is an analysis tool for considering corneal haze resulted from laser eye surgery and is detailed in this section to provide background knowledge of corneal haze. Light scattering takes place in a case that an electromagnetic (EM) wave, the incident light, encounters scattering particles. At the moment the EM waves collide with discrete particles, electrons oscillate within the particles at the same frequency as the incident wave. The oscillation, called a dipole moment, is regarded as a source of light scattering. The energy of the incident light is either discharged by light radiation or extinguished by absorption within the particle. The frequency of the incident light is PAGE 27 14 equal to that of scattered light for what is considered elastic scattering, while Raman scattering is considered an inelastic scattering process. Figure 1-3 shows the light scattering response to an incident light. Electromagnetic Wave + Scattered Light Figure 1-3. Light scattering response to an incident light. There are two kinds of categories for the elastic light scattering. One is the theory of Rayleigh scattering that is applied to a system with small, dielectric (non-absorbing)and spherical particles. The other is the theory of Mie scattering that is used for general spherical solution without a particular size limit. Hence, Mie scattering has no size limit other than geometrical optics so that it can be used for describing most spherical scatterer systems, including Rayleigh scattering. However, Rayleigh scattering is usually used as long as it is applicable due to complexity of Mie scattering. The criteria of Rayleigh scattering is that <<1 and m<<1, where is the dimensionless size parameter given by 2a , (1-1) where a indicates the particle radius, and is the relative wavelength defined as PAGE 28 15 oom , (1-2) where o means the incident wavelength in vacuum and mo represents the refractive index of the surrounding medium. The refractive index, the property of the material, is defined as mni . (1-3) In this equation, n indicates the common refraction of light while the complex term relates to absorption. The value of k is not exactly zero for any material, but materials with the value approaching to zero are termed dielectrics. The relative refractive index is defined as ommm . (1-4) The magnitude of the relative refractive index, m is 1222onmm . (1-5) In this study, Mie scattering theory was applied and described in accordance with the treatment and notation of Kerker (1969). Furthermore, the scattering of light characterized by the Mie scattering solution for single particles was fulfilled by following three criteria reported by Jones (1979). The scattering particles must have enough space not to occur significant electrical interaction between them, which is satisfied if the distance between scattering particle centers is two to three times greater than the diameter of particles (Kerker, 1969). Only the smallest mean free pathlength for the corneal simulation, 5 m, may conflict with this PAGE 29 16 requirement. However, the 5 m pathlength is probably not clinically relevant (see Table 2-2). Figure 1-4 shows the geometry of a Mie light scattering particle relevant to the following equations. Figure 1-4. Mie light scattering particle geometry. For each scattering angle (,), the equations 1-6 and 1-7 represent the intensities of scattered radiation vertically and horizontally polarized with respect to the scattering plane 22122sin4oIIir , (1-6) 22222cos4oIIir , (1-7) where Io is the incident intensity ,and the intensity functions i1 and i2 are defined as 21121coscos1nnnnnniabnn , (1-8) PAGE 30 17 22121coscos1nnnnnniabnn . (1-9) In the equations 1-8 and 1-9, the angular dependent functions n and n Legendre polynomials are (1)coscossinnnP , (1-10) (1)coscosnndPd , (1-11) and where the parameters an and bn are defined as '''nnnnnnnmmmammm 'n , (1-12) ''''nnnnnnnnmmmbmmm . (1-13) The size parameter is 2ooam . (1-14) The equation 1-15 represents the parameter n, the half-order Bessel functions of the first kind (Jn+1/2(z)). 12122nnzzJ z . (1-15) The equation 1-16 describes the parameter n 12122nnnzzHzzi nz , (1-16) PAGE 31 18 where H(2)n+1/2(z) is the half-order Hankel functions of the second kind. The parameter Xn is the half-order Bessel functions of the second kind Yn+1/2(z), which is given by 12122nzzY nz . (1-17) The vertical-vertical differential scattering cross section indicated in the equation 1-18 means that both the incident light and the scattered light after the interaction are vertically polarized with respect to the same scattering plane (yz-plane) 2'124VVC i . (1-18) Similarly, the horizontal-horizontal differential scattering cross section shown in the equation 1-19 means that both the incident light and the scattered light are polarized parallel to the scattering plane (xz-plane) 2'224HHC i . (1-19) In the equations 1-18 and 1-19, the first subscript means the incident light, and the second subscript means the scattered light. Also, subscripted V and H respectively refer to the vertical polarization and the horizontal polarization. The intensity of the scattered light in any direction for the unpolarized incident light is given by 212224o I Iir i , (1-20) and, the unpolarized differential scattering cross section is given by 2'1228scatCi i . (1-21) The unit of cm2/sr is used in the present study for differential scattering cross-sections. PAGE 32 19 1.4 Monte Carlo Simulations Monte Carlo (MC) simulation is a numerical technique that could statistically solve complex problems by means of a computer simulation. It has been applied to a variety of research fields since it was first proposed in 1949. It utilizes probability information on any number of parameters to approach a problem that would be hard or impossible to be solved with other methods (Niemz, 1996). Above all, the main advantage of this simulation is that the simple and versatile application is possible. It enables specific boundary conditions and source/detector locations easily to be added to wherever possible. MC simulations have been used to study light scattering in tissue and to grasp light tissue optical properties, which is essential in biological laser applications. In a study of a breast tissue imaging, a model of multiple light scattering was simplified in accordance with the linear transport theory. The Monte Carlo sampling was controlled by relevant probability distributions and photon transport theory, and individual photons were tracked as they moved through simulated breast tissue. It was turned out that a problem of this simulation is the trade-off between precision and computing time (Key et al., 1991). In another study, MC simulations were successfully used to verify experimental results of skin perfusion conducted by a laser Doppler blood flowmeter. The simulations were also used for analyzing different wavelengths and optical probe geometries, and a system which could probe skin at varying depths was developed in this study (Koelink et al., 1994). PAGE 33 20 In a polystyrene spheres suspension study, the MC simulation was used to confirm the scattering coefficient and anisotropy factor of blood obtained by an experiment with laser Doppler measurements. (Kienle et al., 1996). In a study on the human brain, MC simulations were also used to determine optical properties, including scattering and absorption coefficients (Bevilacqua et al., 1999). MC simulations have been used to develop a laser induced breakdown spectroscopy for continuous emissions monitoring (Hahn et al., 1997). MC simulations were also used to solve the linear Boltzmann equation of photon transport for varying mean-free paths, and then the results from the simulation were compared with experimental results of the polarization-dependent scattering matrix using polystyrene latex spheres (Kaplan et al., 2001). MC simulation facilitating a better understanding of the photo physical interaction of excitation and fluorescent light in tissue was used by Pogue and Burke to design a tissue fluorescence probe and to confirm their experimental results (Pogue & Burke, 1998). 1.5 Objectives The purpose of this study is to develop the objective method that is a two-color ratio technique to quantify corneal haze clinically. Previously, the two-color ratio technique study was conducted using PS suspensions coupled with animal studies with eyes of bovine and rabbit, and the MC simulation (Coffey, 2001). In this study, the experimental systems and methods were improved to better couple with the optimal simulation, and the MC simulations were extended to a three dimensions. Specifically, this study focuses on improving the outcome and methodology of the previous study. PAGE 34 21 This study is primarily aimed to develop the three-dimensional Monte Carlo computer program; therefore, the cornea could be simulated more ideally, and the simulation can yield more reasonable outcome. Based on the computer simulation, the fiber optic probe will be built so that it could contribute to obtaining data at desired positions in the PS suspensions experiments. With the results through both PS suspensions experiment and MC simulation, the two-color scattering ratio will be calculated and verified. Ultimately, the promising two-color light scattering technique will be further developed in support of the goal of the in vivo measurement of the haze and a consistent quantitative metric for the haze assessment of the human cornea. Therefore, it will help research in the healing of corneal tissue and haze formation. PAGE 35 CHAPTER 2 POLYSTYRENE SUSPENSION EXPERIMENTS For the purpose of comprehending scattering phenomena, experiments were conducted using polystyrene (PS) suspensions with known scattering particle density and size. In the experiments, the PS suspensions were used to simulate the dense scattering environment of the human cornea connected with clinical haze. 2.1 The Three-Color Laser System Experiments on the PS suspension were conducted with a three-color laser system. The system was designed to compare the intensity levels of three laser wavelengths at each collection fiber optic. The three lasers employed were a helium: neon laser (wavelength of 594 nm) and two diode lasers (wavelength of a 532 nm and 670 nm). The components of this laser system were mirrors and wedges, including the three lasers. With respect to the overall three-color laser system, a schematic and a picture of the experimental setup are respectively shown in Figures 2-1 and 2-2. The three lights initially emitted from the lasers are transmitted through the optical system, brought together coaxially, and then launched into the illumination fiber optic. The function of the illumination fiber optic is the transmission of the laser light into a PS suspension. Four fiber optics, including an illumination optic, with a diameter of 400 m are fixed through holes in an aluminum bar. The diagram of the fiber optic and the placement of fiber optics in the aluminum bar are shown in Figures 2-3 and 2-4. 22 PAGE 36 23 Figure 2-1. Picture of experimental setup of the three-color laser system. PAGE 37 24 Figure 2-2. Schematic of experimental setup of the three-color laser system. Figure 2-3. Diagram of the fiber optic. PAGE 38 25 560um 400um 160um 200um Figure 2-4. The placement of fiber optics in the aluminum bar. Three collection fiber optics were sequentially placed next to the illumination fiber optic at intervals of 560 m so that they could collect and transmit the scattered light to a spectrometer/detector directly attached to them. Scattering signals of three laser lights were measured synchronously with an integrated spectrometer/detector unit. The detector was a 2048-element linear CCD array with an optical resolution of approximately 0.69 nm and effective dispersion of 0.13 nm/pixel. The spectrometer/detector was controlled by software, with integration time of 100 ms for this study. Apparatuses for the three-color laser system are listed in Table 2-1. 2.2 Experimental Method The similar experimental procedures and PS suspensions which were made in the previous study (Coffey, 2001) were used in this study. Moreover, preliminary testing for eliminating the external environment effects and variations performed by Coffey were accepted. These were summarized below. PAGE 39 26 Table 2-1. List of apparatuses for the three-laser system. Device Manufacturer Model Description Laser and Electronics 532 nm Diode Laser Lasermate Group, Inc GME-532-3XF9P-Z1 532 nm, 1.5~3 mW 594 nm Helium:Neon Laser JDS Uniphase 1677 594 nm, 4 mW 670 nm Diode Laser Edmund Industrial Optics LDL 175 670 nm, 5 mW Fiber Optic Spectrometer Ocean Optics, Inc S2000 Production 500-770 nm, 1200 line/mm, 10-m Silt Laptop Computer Gateway SOLO 5300 Intel Celeron Ocean Optics OOlbase32 Optics Fiber Optic Thor Labs, Inc FG-200 LCR 1 Meter Mirror 1 Newport Corp. 10D10ER.1 25.4 mm Flat Pyrex Mirror Mirror 2 Newport Corp. BB 1-E02 Dielectric Mirror Parallel Windows Newport Corp. 10QW20-30 25.4 mm Uncoated Glass Wedge 2.2.1 Polystyrene Suspensions Clinically pertinent PS suspensions were created based on information from the literature which was used to estimate the scattering particle (damaged keratocytes) size and the number density of the particles at varying times during the healing process after laser refractive surgery (Fantes et al., 1990). Typically, haze becomes clinically relevant about seven days post-operative and is maximized around three weeks post-operative. In Table 2-2, relevant scattering parameters were summarized referring to the research reported by Fantes and co-workers concerning wound healing in monkeys. Table 2-2. Scattering particle density and mean scattering pathlength at 7-days and 3-weeks after laser ablation. Scattering Density (Fibroblasts/ml) Scattering Pathlength (m) 7 days 53,000,000-76,000,000 29-33 3 weeks 303,000,000 19 (Fantes et al., 1990) PAGE 40 27 The corresponding light scattering pathlength values (i.e. mean free path) between the scattering particles were calculated by Equations 2-1 and 2-2, 163xNfv , (2-1) where, lx , (2-2) fv the total particle volume fraction N the particle number density x the diameter of the space surrounding each scattering particle l the mean scattering pathlength. Each scattering particle is surrounded by a spherical region, and the overall volume fraction is assumed as a unity in order to simulate equal distance of separation among scattering particles as shown in Figure 2-5. In the case of light scattering, the diameter of the space surrounding each particle(x) is equal to the pathlength (l). Scattering Particle d x l Figure 2-5. Expression of scattering particles in suspension as a volume fraction of one. The suspensions were created with monodisperse PS microspheres of 7.20 m in diameter from Bangs Laboratories, Inc. The microspheres were liquid suspensions with PAGE 41 28 10.0% solid content by weight and a density of 1.062 g/ml. The specifications of the microspheres are detailed in Table 2-3 for the 7.20 m particle. Table 2-3. Polystyrene microsphere specifications. Particle Size 7.2 m Product Code PS06N/1343 Standard Deviation (m) 0.07 Microspheres/ml 4.846E+08 Based on the clinical haze conditions and serial dilution, the original stock suspensions were diluted with ultra-purified deionized (DI) water to create four different concentrations which were kept in 15 ml glass bottles and named as C1, C2, C3 and C4. The densities of the PS suspensions are specified in Table 2-4, as based on the dilution factors. Table 2-4. Polystyrene suspensions densities for particle size of 7.20 m. PS Suspension Density (Microspheres/ml) Volume Fraction (fv) (%) C1 96,920,000 1.27 C2 48,460,000 0.63 C3 9,692,000 0.13 C4 969,000 0.01 The following four suspensions were prepared for the particle size of 7.20m, and were stored in the refrigerator. C1: the most dense suspension C2: the second most dense suspension C3: the third most dense suspension C4: the least dense suspension In these suspensions, the particles were observed to settle out of suspension and aggregate on the bottom of the bottle over time. All visible signs of aggregation were PAGE 42 29 readily removed through multiple cycles of mechanical vortexing and sonic agitation before measurements. 2.2.2 The verification of Suspension Concentration Transmission measurements were conducted to verify the actual concentration of each of the four suspensions. Figure 2-6 shows a schematic of experimental setup of transmission measurement. Figure 2-6. Schematic of transmission measurement. A 632.8 nm helium: neon laser was used in this measurement, and the aperture was used for blocking forward scattered light. Reference measurements were performed with the purified DI water and the intensity (I0) was measured with power meter. Similarly, the intensity of each suspension (I) was measured with power meter. Measurements of each PAGE 43 30 suspension were performed twelve times, and reference measurements were recorded before and after suspension measurements. As denoted in the following equation 2-3, the transmission () was determined by dividing the intensity of each suspension by the average of the reference intensity (i.e. DI water only). afterbeforeIII002 . (2-3) The mean and standard deviation values of transmission for each suspension are enumerated in Table 2-5. Table 2-5. Mean and standard deviation values of transmission for each suspension. PS Suspension Mean Standard Deviation C1 0.033102 0.001291227 C2 0.234376 0.003564348 C3 0.360230 0.007040336 C4 0.312850 0.001285893 It was noted that the concentration of some suspensions was too high to measure the transmission with the power meter; therefore, the suspensions were diluted by adding the DI water. The suspensions of C1 and C2 were diluted by adding suspension of 0.1 ml to the DI water of 2.9 ml. For the suspension C3, 1 ml of suspension was added to 2 ml of DI water in six measurements, and 0.2 ml of suspension was added to 2.8 ml DI water in the other six measurements. On the other hand, suspension C4 was not diluted. With calculated transmission and dilution factor (f), the extinction coefficient () was calculated by the equation below, where L is the 1 cm pathlength, extK LfKext)ln( . (2-4) Equations 2-5 through 2-7 express how the standard deviation of the extinction coefficient was determined. The standard deviation of the extinction coefficient is defined PAGE 44 31 in Equation 2-5, where the derivatives are expressed in Equations 2-6 and 2-7. In these equations, represents the standard deviation of the transmission, and f represents the standard deviation of the dilution which is assumed as 5% of the mean dilution. 22ffKKextextKext (2-5) fKext1 (2-6) 2lnffKext (2-7) Table 2-6 lists values of the dilution factor, the extinction coefficient and the relative standard deviation for each suspension. Table 2-6. Values of the dilution factor, the extinction coefficient and the relative standard deviation for each suspension. Standard Suspension Dilution Kext Deviation Factor (f) of Kext C1 0.0402 88.8232 4.3512 C2 0.0333 43.6008 2.2236 C3 0.2000 8.6910 0.2733 C4 1.0000 1.0229 0.0582 The actual number density (N (particles/cm3)) in suspensions was calculated by the equation 2-8, LCKNextext , (2-8) where Cext is the extinction cross-section, which denotes total abstraction of incident light from the scattering particle because of both scattering and absorption. Using the DBMIE code, the total extinction cross-section, 9.33e-7 cm2, was calculated for the PS PAGE 45 32 suspensions, using a particle diameter of 7.2 m, refractive index of m=1.2-0.00i, medium index of m0=1.33, and wavelength of 632.8 nm. The mean scattering pathlength of each PS suspension was calculated by using the equation 2-1 and 2-2. Values of number densities and mean scattering pathlength are listed in Table 2-7. It is noted that the number densities are in excellent agreement with the Table 2-4 values, which are based on serial dilution of the stock suspension. Table 2-7. Values of number densities and mean scattering pathlength of suspensions. PS Suspension Number Densities Mean Scattering Pathlength C1 9.5e7 27.2 C2 4.7e7 34.4 C3 9.3e6 59.0 C4 1.1e6 120.3 2.2.3 Measurement Method In order to control variations in laser power, reference measurements were performed with a small block of commercial-grade ultra-high molecular weight polyethylene (PE) before and after PS suspension measurements. All measurements were carried out perpendicular to the surface by holding the fiber optics in direct contact. For a fixed angle and a precise depth in the PS suspensions, fiber optics were held by a vertically mounted micrometer stage while taking measurements. The PS sample bottles were held in place by a block of aluminum with a hole machined to fit the sample bottle diameter. Figure 2-7 represents pictures of the measurements for PE reference and PS suspension. Figure 2-8 displays a representative scattering spectrum of PE measurement. PAGE 46 33 Figure 2-7. Pictures of measurements for PE reference (Left) and PS suspension (right). 020040060080010001200500550600650700750IntensityWavelength (nm)670-nm Peak594-nm Peak532-nm PeakBaseline Figure 2-8. A representative scattering spectrum of PE measurement. PAGE 47 34 2.2.4 The Scattering Sample to PE Reference Ratio The experimental variations because of day-to-day changes of laser intensity would yield unwanted results. Such experimental problems can be eliminated in the signal-to-reference ratio that was calculated with the absolute peak intensity of spectra. The absolute peak intensity of each recorded spectra for each scattering wavelength was initially calculated by subtracting the average baseline intensity from the spectrum intensity. As shown in Figure 2-9, an area marked with oblique lines indicates the absolute peak intensity which was calculated in the equation 2-9, basePeakxxiipeakyWyARIGHTPEAKLEFTPEAK . (2-9) In this equation, W and Peak basey are respectively defined as the equations 2-10 and 2-11, 1 LEFTPIXELRIGHTPIXELPeakxxW , (2-10) 2)()(RightbaseLeftbasebaseyyy . (2-11) The scattering sample-to-PE reference ratio was then calculated by dividing the absolute peak area from the PS scattering sample by the average of two PE reference peaks measured before and after the sample measurement. The equation 2-12 defines the scattering signal-to-reference ratio ( R ), )(21afterbeforePEPEPSR , (2-12) where, PS the PS spectra absolute scattering peak area, PEbefore and PEafter the PE spectra absolute scattering peaks. PAGE 48 35 Figure 2-9. The absolute peak intensity. 2.2.5 Preliminary Testing It was necessary to check whether effects from the external environment of the PS suspension would have an influence on the experimental results. For external effects attributable to the glass vials and laser intensity variations, preliminary testing was performed in the previous study (Coffey, 2001). Briefly, it turned out that the light reflected from the glass vials did not cause a significant influence on the experiments results due to the depth of the vials. Also, laser intensity variations and integration time effects are canceled out in the scattering signal-to-reference ratio as discussed before. 2.3 Data Analysis The PS suspensions were analyzed to determine the light scattering signal (signal-to-reference ratio) dependence on particle concentration and scattering wavelength. To do this, scattering measurements were recorded for each PS suspension. All spectra were PAGE 49 36 recorded from 12 averages (recorded on different days) of individual 100 ms integration times. For each experiment, spectra of PS suspension were recorded, and PE spectra were recorded before and after recording PS spectra. All measurements were taken 0.5 mm below the suspension surface. The resulting signal-to-reference values for each wavelength and concentrations at each optic, along with their standard deviations are tabulated in Table 2-8 and shown in Figure 2-10 through 2-12. Table 2-8. The PS signal-to-PE reference ratios for each wavelength and concentrations at each fiber optic. Sigma represents the standard deviation of the signal mean (N=12). Fiber optic 2 Fiber optic 3 Fiber optic 4 Wavelength Suspension Pathlength PS/PE Sigma PS/PE Sigma PS/PE Sigma C1 27.2 1.77 0.40 1.24 0.32 1.13 0.33 532 nm C2 34.4 1.04 0.18 0.64 0.11 0.53 0.11 C3 59.0 0.53 0.03 0.22 0.02 0.18 0.01 C4 120.3 0.44 0.18 0.15 0.02 0.12 0.01 C1 27.2 2.18 0.44 1.58 0.39 1.35 0.40 594 nm C2 34.4 1.33 0.19 0.84 0.16 0.67 0.14 C3 59.0 0.64 0.04 0.28 0.02 0.22 0.01 C4 120.3 0.56 0.24 0.19 0.02 0.15 0.02 C1 27.2 3.50 0.67 2.41 0.58 2.05 0.56 670 nm C2 34.4 2.21 0.35 1.36 0.24 1.07 0.23 C3 59.0 1.18 0.15 0.53 0.06 0.41 0.04 C4 120.3 1.04 0.48 0.37 0.06 0.29 0.07 For all cases, the scattering signal-to-reference ratio increases as the concentration of suspensions increases. The greater the quantity of scattering particles in the suspension is, the greater the quantity of light is scattered back into the collection optics. PAGE 50 37 0.000.501.001.502.002.503.003.504.004.500.020.040.060.080.0100.0120.0140.0ConcentrationSignal-to-Reference Ratio 532nm 594nm 670nm Figure 2-10. The scattering signal-to-reference ratios at fiber optic 2. The error bars denote standard deviations12 sets of experimental mean performed on different days. 0.000.501.001.502.002.503.003.500.020.040.060.080.0100.0120.0140.0PathlengthSignal-to-Reference Ratio 532nm 594nm 670nm Figure 2-11. The scattering signal-to-reference ratios at fiber optic 3. The error bars denote standard deviations12 sets of experimental mean performed on different days. PAGE 51 38 0.000.501.001.502.002.503.000.020.040.060.080.0100.0120.0140.0PathlengthSignal-to-Reference Ratio 532 nm 594 nm 670 nm Figure 2-12. The scattering signal-to-reference ratios at fiber optic 4. The error bars denote standard deviations12 sets of experimental mean performed on different days. The two-color scattering ratios were also calculated from the concentration analysis data. The color ratio, CR, is defined in Equation 2-13. It is the signal-to-reference ratio of a lower wavelength divided by that of a higher wavelength, 21RRCR . (2-13) The resulting two-color ratio values for the laser wavelength combinations of 532 nm/594 nm, 594 nm/670 nm, and 532 nm/670 nm and PS concentrations at each fiber optic, along with their respective standard deviations are tabulated in Table 2-9 and shown in Figures 2-13 through 2-15. PAGE 52 39 Table 2-9. Two-color ratios for each concentration at each fiber optic. Sigma represents the standard deviation of the signal mean. Fiber optic 2 Fiber optic 3 Fiber optic 4 Wavelength Suspension Pathlength Ratio Sigma Ratio Sigma Ratio Sigma C1 27.2 0.81 0.05 0.79 0.02 0.83 0.09 532/594 C2 34.4 0.79 0.05 0.77 0.03 0.79 0.02 C3 59 0.82 0.05 0.79 0.03 0.81 0.03 C4 120.3 0.81 0.13 0.79 0.06 0.77 0.07 C1 27.2 0.62 0.05 0.65 0.04 0.66 0.09 594/670 C2 34.4 0.6 0.07 0.62 0.05 0.62 0.03 C3 59 0.55 0.04 0.54 0.04 0.55 0.05 C4 120.3 0.52 0.03 0.51 0.04 0.52 0.05 C1 27.2 0.5 0.06 0.51 0.04 0.55 0.05 532/670 C2 34.4 0.47 0.07 0.47 0.04 0.49 0.03 C3 59 0.45 0.05 0.42 0.03 0.45 0.04 C4 120.3 0.42 0.09 0.4 0.05 0.4 0.07 0.300.400.500.600.700.800.901.000.030.060.090.0120.0150.0PathlengthTwo Color Ratio 532/594 594/670 532/670 Figure 2-13. The two color ratios at fiber optic 2. The error bars denote standard deviations12 sets of experimental mean performed on different days. PAGE 53 40 0.300.400.500.600.700.800.900.030.060.090.0120.0150.0PathlengthTwo Color Ratio 532/594 594/670 532/670 Figure 2-14. The two color ratios at fiber optic 3. The error bars denote standard deviations12 sets of experimental mean performed on different days. 0.300.400.500.600.700.800.901.000.030.060.090.0120.0150.0PathlengthTwo Color Ratio 532/594 594/670 532/670 Figure 2-15. The two color ratios at fiber optic 4. The error bars denote standard deviations12 sets of experimental mean performed on different days. PAGE 54 41 The color ratio increases with increased scattering concentration (increased light scatter), with the single exception of the 532 nm/594 nm ratio. It can be concluded that there is a unique color dependence of light scattering under these dense scattering medium conditions because the color ratio changes with the concentration of scattering particles. However, the calculated absolute values of two-color ratio cannot be considered as the pure values. As shown in equation 2-14, the calculated values still include the PE reference values that were measured to eliminate the variation due to laser power change. In this equation, 21CC represents the two-color ratio for PS suspension, while '2'1CC represents the two-color ratios for the PE reference material. The PE reference values must be eliminated from the two-color ratios to obtain the pure values, which is discussed in Appendix C. '1'221expCCCCR . (2-14) PAGE 55 CHAPTER 3 MONTE CARLO SIMULATIONS Monte Carlo simulations were used to provide solutions for the complex problem of light scattering from abnormal keratocytes in a dense tissue medium and to give support to the experimental works by performing numerical light scattering simulations on a computer. The Monte Carlo computer program tracks an individual photon traveling in a simulated dense scattering medium such as the human cornea. A photon, initially launched into the medium undergoes two types of interactions, which are either absorption or scattering. That is, an individual photon may be absorbed by the fibroblastic keratocytes in the cornea, or it is scattered from the fibroblastic keratocytes until it escapes from the cornea. 3.1 Development of Overall Algorithm The overall algorithm is made up of launching a photon according to an initial launch angle and a rotational angle, and then tracking the photon that scatters according to Mie light scattering theory and progresses in selected distances between scattering centers. Consideration of the rotational angle was newly added to the previous 2-D study (Coffey, 2001) to achieve a new three-dimensional simulation; therefore, four probability sample arrays named the scattering angle (), the launch angle (), the rotational angle () and the pathlength (x) were created in the algorithm. Moreover, it was attempted to make the current simulation more ideal by adding the condition of an absorption effect. The explanation of the development of the probability distributions for these four 42 PAGE 56 43 parameters follows below and a complete copy of the FORTRAN code is attached in Appendix B. 3.1.1 Probability Distribution of Scattering Angle The probability distribution of the scattering angle was defined based on Mie light scattering theory. A photon scatters from a scattering center as shown in Figure 3-1. Negative Direction Incident Light Scattering Particle Positive Direction Figure 3-1. Diagram of a light scattering from a scattering particle. In order to set up a probability array for the scattering angle theta, the first step is to calculate both the vertically polarized differential scattering cross section (Câ€™VV) and the horizontally polarized differential scattering cross section (Câ€™HH) for each integer scattering angle from 0 to 180 degrees. The FORTRAN code DBMie (DAVE, 1968) was used to do these calculations with input values such as the scattering particle size, relative wavelength and relative refractive index of the scattering particles. The relative refractive index was given as 1.08 â€“ 0.01i for representative of cornea. The real part of the refractive index was referred to a study (del Val et al., 2001) while image part was approximated to represent some weak absorption, in that keratocytes are not believed to PAGE 57 44 be true dielectrics. Moreover, a scattering particle diameter (d) was given 6.30 m for varying wavelengths: 532 nm, 594 nm and 670 nm. These conditions represent typical characteristics of light scattering in a human cornea involved with clinical haze. The relative wavelength and relative refractive index are determined as ,oom (3-1) ,ommm (3-2) where m0 and 0 represent the index of the scattering medium and the scattering wavelength in vacuum, respectively. According to Mie light scattering theory, the unpolarized differential scattering cross section (Câ€™UP or Câ€™scat) was defined as )(21'''HHVVscatCCC . (3-3) Figure 3-2 shows the representative differential scattering cross section as a function of angle and wavelength for m =1.08 â€“ 0.01i and d=6.3 m for three different wavelengths. Note the resonant structure characteristic of Mie scattering and the nearly 7-order of magnitude variation in the differential cross section. As shown in Figure 3-2, the scattering response is dependent on wavelength, and it is also noted that forward scattering for Mie-sized (d>>) scattering particles is prominent. With the scattering angle distribution array created on a basis of Mie light scattering theory, a sampling array was built for given values of scattering particle size, relative refractive index and relative wavelengths. The differential scattering cross PAGE 58 45 1101001000100001000001000000100000001000000000306090120150180Scattering Angle (Degree)Differential Scattering Cross Section 532nm 594nm 670nm Figure 3-2. Differential scattering cross section for three wavelengths. sections calculated were very small values ranging from 10e-5 to 10e-11; therefore, it was necessary to normalize those values. As the first step of the normalization, the total differential scattering cross section was calculated by summing all differential scattering cross sections in the array. Second, an array of the normalized differential scattering cross sections was created with each differential scattering cross section divided by the total differential scattering cross section value. Third, the normalized differential scattering cross sections were multiplied by the inverse of the minimum value of the normalized differential scattering cross sections in the array; thereby, the minimum value of the normalized differential scattering cross sections of the array becomes one. After the process of the normalization, all values in the array were converted into integer values, and a sampling array was generated by two kinds of approaches. The overall process of creating the sampling array is indicated in Figure 3-3. PAGE 59 46 Set values in theta (THETD(I)) and degree(DEG(I)) to 0 to 180 Calculate CVV and CHH values for angles 0 to 180 Calculate Cunpolarized (Cup) and the sum of Cup (CTOT) Normalize the Cup array to Cnormalized value (CNORM) CNORM=Cup/CTOT Create the scaled array CNM (Cnormalized and multiplied) CMIN=MINVAL (CNORM) CNM=CNORM/CMIN Convert CNM into an integer array (NM) NM (I)=ANINT(CNM(I)) The Previous Approach The New Approach Build the probability array(PROBAR) for angle selection PROBAR (K)=I-1 Build a new array(PPSUM(I)) for sum of an Integer array(NM(I)) PPSUM (I)=PPSUM(I)+NM(I) Random sampling of THETA SPOT=ANINT(ITOTAL*(1.0-RN)) IF (PPSUM (I) PAGE 60 47 In the approach that was used in the previous study by Coffey, the sampling array was created by sequentially placing each scattering angle from 0 to 180 into the array the number of times equivalent to its corresponding normalized cross section value. For instance, each angle of 179 and 180 was placed into the sampling array one time and two times respectively in a case that the normalized values of the scattering angle distribution array are one and two each for its corresponding angles. The sampling array was accessed at random in order to randomly define scattering angles in proportion to the angular scattering distribution according to the overall Mie theory. Meanwhile, a new approach was attempted in this study for the purpose of optimizing the program run time, which was one of the future works suggested in the previous study. Instead of building a large scale sampling array, a small scale array was built by sequentially placing cumulative values of each of the normalized cross sections into its elements corresponding to the scattering angles, rather than sequentially constructing it the number of times equal to each normalized cross section value for the scattering angles. For example, in case that the values of the normalized cross section were ten and five respectively for the angle of 0 and 1 degree, ten and fifteen would be placed into the first and second element of the sampling array in a row. Consequently, modifying the code enabled to save the code run time by an average 14 seconds for 100,000,000 in the trial, which provided a great optimization of run time for the whole simulation. Furthermore, the new code did not require a slight absorption factor (k of 0.01) in the refractive index that was necessary to create the sampling array because of FORTRAN restrictions on array size in the previous code. Without slight absorption, the PAGE 61 48 resonant structure of the Mie cross-secti ons yielded a longer normalization factor, and thereby generated errors by exceeding the allowable array size limits. 3.1.2 Probability Distribution of Launch Angle The probability distribution of the launch angle was created based on a Gaussian distribution. The angular intensity of the laser light exiting the fiber optic probe tip was assumed to be normally distributed. The distribution of the laser light emitted from the fiber optic is shown in Figure 3-4. The mean angle was 0.0 , and a standard deviation of 2.33 was calculated in accordance with Empirical rule ( 3 ). Generally, the Gaussian probability distribution is defined as 2 21 exp, 2 2 p (3-4) where the expression of the geometric mean ( ) is shown in Equation 3-5, 2 2 p d . (3-5) Figure 3-4. Distribution of the laser light exiting from the fiber optic probe tip. PAGE 62 49 The center of the laser light has the strongest intensity whereas the intensity is the weakest at the edge of the laser light; therefore, the maximum probability of is the same as that of mean. From equation 3-4, the probability of mean is 1.2meanp (3-6) The launch angle probability distribution array was built on a condition that a probability of the minimum (min) and maximum (max) angles of exiting laser light is equal to 1% probability of mean (mean) angle. Equation 3-7 was used to calculate the 1% probability, 2min,maxmaxmin2111exp.100222pp (3-7) Equations 3-8 and 3-9 show the expressions for the maximum and minimum alpha, 01.0ln22max , (3-8) 2min2ln0.01 . (3-9) The launch angle sampling array was then created in accordance with the launch angle probability distribution array. The alpha ranged from min to max was divided into 100 increments, and the probability at each increment was calculated. For each increment 1 to 101, the integer value of alpha was placed into the sampling array the number of times equivalent to its percentage probability. 3.1.3 Probability Distribution of Rotational Angle The probability distribution of rotational angle was generated based on uniform distribution. As shown in Figure 3-5, the rotational angle is how much the scattering plane is rotated about z-axis. The angle is uniformly distributed from 0 to 360 in the PAGE 63 50 perpendicular plane with respect to the scattering plane. To describe a photon traveling in the three-dimensions, a new position of a photon initially scattered or launched in the two dimensional scattering plane was determined by taking a rotational angle into consideration. The probability of rotational angle is the same for each angle from 0 to 360 degree; therefore, the mean angle and standard deviation are all 0.0 . The probability distribution array was created on the condition that all angles had equivalent probability. Figure 3-5. Representation of the rotational angle. The rotational angle sampling array was then created in accordance with the rotational angle probability distribution array. The range of beta from min to max was divided into 360 increments, and the probability at each increment was calculated. For each integer 1 to 361, the incremental value of beta was entered into the sample array the number of times equal to its percentage probability. x z y PAGE 64 51 3.1.4 Probability Distribution of Pathlength The probability distribution of pathlength was produced based on a lognormal distribution, which does not allow the negative values of pathlength between scattering particles. The mean pathlength was chosen an input parameter, and the standard deviation was calculated by dividing the value of the mean pathlength by 3.75 on the assumption of observations in corneal histology. Generally, the lognormal probability distribution (Espenscheid, 1964) is defined as 2222lnlnexp22omoomxxepxx . (3-10) The mean pathlength ( x ) and standard deviation () of the pathlength are shown in equation 3-11 and 3-12, which are defined in terms of the modal pathlength () and a dimensionless measure of variance (0) of the pathlength. mx 322lnln,mxx o (3-11) 2exp4exp3.mox 2o (3-12) Equations 3-13 and 3-14 express 0 and of the pathlength obtained from Equations 3-11 and 3-12. mx 1ln2xo , (3-13) 223expomxx . (3-14) PAGE 65 52 The maximum probability is the same as the probability of the modal value (); therefore, the maximum probability is mx 22.2omomepxx (3-15) The pathlength probability distribution array was built by using 1% of the maximum probability. Equation 3-16 was used for calculating 1% probability to create the minimum () and maximum () values in the array, minx maxx 22222min,maxmaxmin2lnln1exp.100222oomoomomxxeepxpxxx (3-16) Equation 3-17 and 3-18 represent the minimum and maximum pathlength, 01.0ln2exp2maxomxx , (3-17) 01.0ln2exp2maxomxx . (3-18) The pathlength sampling array was then created in accordance with the pathlength probability distribution array. The range from xmin to xmax was divided into 100 increments, and the probability at each increment was calculated. For each integer 1 to 101, the incremental value of pathlength was placed into the sampling array the number of times equivalent to its percentage probability. 3.1.5 Absorption The individual photon traveling in the medium goes through two types of interactions, namely scattering or absorption. A plot over the transmission of human cornea defined by Sacks at al (1998) was used to determine the absorption coefficient. As PAGE 66 53 shown in Figure 3-6, the transmission is approximately 92% for three wavelengths, 532 nm, 594 nm and 670 nm, respectively. The general equation for the transmission and absorption coefficient is )exp(0LKIIabsL , (3-19) Figure 3-6. Transmission of the human cornea. where is the transmission, and is the absorption coefficient. The length (L) from anterior surface to posterior surface of the cornea was assumed 500 m constantly. The absorption coefficient can be obtained through the equation 3-19; however, the transmission of 92% is measured in vitro. Because the transmission in vivo is needed, the different factor () between in vivo and in vitro was calculated by equation 3-20. Also, the transmission of 0.958 was given by Equation 3-21; hence finally an absorption coefficient of 8.595E-5 (1/m) was calculated by Equation 3-19. absK 33.1;12121212nnnnnn , (3-20) PAGE 67 54 2)1(92.0 . (3-21) Before a new position of the individual photon was determined, the transmission was calculated with the absorption coefficient and scattering pathlength randomly selected. In case that transmission is greater than a certain scale determined by the Random Number Generator, a photon keeps traveling until the next interaction occurred. On the other hand, it is assumed that a photon is absorbed by a scattering particle in case that transmission is less than the scale. As mentioned above, the transmission calculated with 500 m of pathlength is 0.958, and the transmission even increases as scattering pathlength decrease. Therefore, the probability that the RNG will be determined a scale larger than the transmission is below 5%. In case that mean scattering pathlength is small, the overall interaction in the cornea increase, thereby the probability of absorption would increase. Around 5% of photons were absorbed in the corneal simulation as listed in Table 3-4 (in the section 3.3.1) 3.1.6 Photon Progression The progress of a photon traveling through the given scattering medium is tracked with respect to proper boundaries. Figure 3-7 shows a simplified two-dimensional model of the cornea for the simulation. A photon is initially launched from the fiber optic on the surface of the anterior cornea. Then, the photon is tracked until it disappears in the described medium. That is, tracking of the photon terminates when the photon vanishes because of being absorbed or leaves the cornea through one of three surfaces such as the posterior surface (rinner), the anterior surface (router), or either side surface. PAGE 68 55 Fiber Optic Probe Anterior Coneal Surface 500 m Posterior Coneal Surface Router=8000 m Rinner=7500 m 40.03 Center of Coneal Curvature Figure 3-7. Simple geometry of the simulated human cornea. 3.1.6.1 Coordinate transformation In a three-dimensional coordinate system, coordinates of x, y, and z are calculated by the following equations, sincosLx , (3-22) sinsinLy , (3-23) cosLz , (3-24) where L, and respectively represent pathlength , the rotational angle and the scattering angle. The calculated values of x, y, and z with the above-mentioned equations are still measured in terms of the relative coordinate system. In order to describe the track of a photon in the medium, the coordinates of a new position of the photon should be determined in terms of the absolute coordinate system. In Figure 3-8 indicating point transformation, the coordinates of point C measured in the coordinate system [O], C is obtained by the equation 3-25, PO PAGE 69 56 P R R P PB C A B O A O B O C . (3-25) In this equation, PB C is the coordinate of point C measured in the B coordinate system, and PO B is the coordinate of point B measured in the O coordinate system, which is calculated by the equation 3-26, P R PA B O A O B . (3-26) RO A and RA B represent 3 3 rotation matrices, which determine the relationship between the coordinate system that is typically defined in terms of rotation about x, y, or z axes. Figure 3-8. Depiction of point transformation. As shown in Figure 3-8, the coordinate transformation can be explained as follows: the A coordinate system is initially aligned with the O coordinate system and then rotated degrees about the z-axis, which refers to the rotational angle. Finally, it is rotated degrees about y-axis, which represents the scattering angle. With regard to rotations about the z-axis, the rotation matrix (,zR ) is defined by the equation 3-27, PAGE 70 57 1000cossin0sincos,zR . (3-27) Moreover, the rotation matrix () about the y-axis is defined by the equation 3-28, ,yR cos0sin010sin0cos,yR . (3-28) By multiplying the equation 3-27 by 3-28, the rotation matrix () between the O and A coordinate system can be written as the equation 3-29, ROA cos0sinsinsincoscossinsincossincoscosROA . (3-29) The rotation matrix between two coordinate systems is consistent with the equation 3-29. Thus, the ith rotation matrix between the original coordinate system and ith coordinate system is defined as the equation 3-30, i. (3-30) RRRRRii12312010 Therefore, the (i+1)th coordinates of a photon measured in the absolute coordinate system are calculated by the following equations, iiiiiiiixzRyRxRx1130112011101 , (3-31) iiiiiiiiyzRyRxRy1230122012101 , (3-32) iiiiiiiizzRyRxRz1330132013101 , (3-33) PAGE 71 58 where, 11,iiyx and 1iz are the relative coordinates of the (i+1)th point measured in the (i+1)th coordinate system, and i represents a scalar component of 1st column by 1st row of the ith rotation matrix. 110R 3.1.6.2 The outline of photon progression The outline of photon progression is briefly explained below. 1. Determine the number of photons for trial and the boundary of the medium. 2. Select a lunch angle randomly from the launch angle sampling array based on the Gaussian distribution. 3. Select a rotational angle randomly from the rotational angle sampling array based on the uniform distribution. 4. Select an initial distance randomly from the pathlength sampling array based on the lognormal distribution. 5. Calculate the initial location of the photon and create the rotation matrix (33) for converting the location to absolute coordinates relative to the external center of the medium. 6. Select a scattering angle randomly from the scattering angle sampling array based on the angular Mie scattering distribution. 7. Select a rotational angle and a pathlength randomly as done at the step 3 and 4. 8. Calculate the transmission that depends on the pathlength randomly selected and compare the transmission with a scale value randomly determined by a random number generator. Only if the transmission is smaller than a scale value, jump to step 14. 9. Calculate the new location of the photon in the relative coordinate system. 10. Calculate the new location of the photon in absolute coordinate system based on the rotation matrix. 11. Calculate the radius of the particle relative to the inner and outer radii and the angle of the radius to the vertical center axis. 12. Verify the location of the photon within the bounds of the medium. 13. Repeat steps 6 through 12 until the photon is out of boundary. PAGE 72 59 14. Repeat steps 2 through 13 for the total trial number of photon. 3.1.7 The Whole Algorithm of the program The schematic flowcharts of the whole program are shown in Figure 3-9 and 3-10. Figure 3-9 describes the sample array generation, and Figure 3-10 explains photon progression through the cornea. The complete FORTRAN code is attached in Appendix B. Sample array generation Input mean launch angle and standard deviation of launch angle Input mean rotational angle and standard deviation of rotational angle Input mean pathlength and standard deviation of pathlength Input scattering particle size, wavelength, and refractive index DBMie used to create differential scattering cross section sample array Gaussian sample array for launch angle was created Uniform sample array for rotationalangle was created Lognormal sample array for pathlength was created Figure 3-9. Flow chart of Monte Carlo program algorithm (Sample array generation). PAGE 73 60 Photon progression through the cornea Input the number of photons to be tracked and the inner and outer radii of the boundary Randomly select a launch angle from the launch angle sampling array Randomly select a rotational angle from the rotational angle sampling array Repeat for input number of photons Randomly select a pathlength from the pathlength sampling array Determine the location of the photon in the cornea Randomly select a scattering angle from the differential scattering cross section sampling array Randomly select a rotational angle and a pathlength from each sampling array If photon is still within corneal bounds Calculate a transmission and determine whether a photon is absorbed Determine the location of the photon in the cornea Progression for a photon ends andthe exit position is recorded Yes If photon is not within Corneal bounds No Figure 3-10. Flow chart of Monte Carlo program algorithm (Photon progression through the cornea). PAGE 74 61 3.2 Verification of Program To make sure of the accuracy of the program operation, each probability distribution array and photon tracking were verified in advance before conducting the Monte Carlo simulation studies. 3.2.1 Verification of the Scattering Angle Distribution The verification of the scattering angle distribution was performed by sampling the scattering angle 100,000,000 times. The input values used for verification of the program were a wavelength of 594 nm, a refractive index of 1.08 â€“ 0.01i, and a scattering particle diameter of 6.30 m. Figure 3-11 shows the normalized Mie distribution series and the differential scattering cross section distribution obtained by running the Monte Carlo simulation code. The distribution of the sample series accurately represents the differential scattering cross-section distribution, which is still comparable with the distribution for the Mie series. 1101001000100001000001000000100000001000000000306090120150180Scattering Angle Theta (Degree)Differential Scattering Cross Section (C'scat) Mie Distribution Sample Distribution Figure 3-11. Verification of the scattering angle theta. PAGE 75 62 Therefore, it can be concluded that the code and the Random Number Generation system accurately create the sampling array for the scattering angle and functions properly without any logical error. 3.2.2 Verification of the Launch Angle Distribution The verification of the launch angle distribution was performed by sampling the launch angle 65,000 times. The trial run of 65,000 samples was performed five times for a total of 325,000 samples. The input values for this verification, the mean angle and standard deviation were respectively 0 and 2.33 considering the geometry of the light exiting from the probe tip (Section 3.1.2). Table 3-1 lists the five trial values of the mean and standard deviation along with the averages of them. The average values of the mean (-0.01) and the standard deviation (2.31) correspond excellently with the input values. Table 3-1. The five trials of the mean and standard deviation for sampling the launch angles, along with the average of them. Input Trial 1 Trial 2 Trial 3 Trial 4 Trial 5 Average Mean ( x ) 0.00 0.00 -0.01 0.00 -0.04 -0.01 -0.01 Standard Deviation x 2.33 2.30 2.31 2.30 2.31 2.31 2.31 Figure 3-12 shows the result of the Gaussian probability function with two sets of the simulation results overlaid onto it. The two scaled sampling arrays show good agreement with the exact Gaussian distribution; therefore, it can be concluded that the code and the RNG system accurately create the sampling array for the launch angle and operates properly without any logical error. PAGE 76 63 02468101214161820-10-8-6-4-2024681Alpha(Degree)Probability (%) 0 Gaussian Distribution Trial 1 Trial 2 Figure 3-12. Verification of sampling algorithm for launch angle alpha. Two sets of MC trials of 65,000 samples with the Gaussian probability function. 3.2.3 Verification of the Rotational Angle Distribution The verification of the rotational angle distribution was performed by sampling the rotational angle 65,000 times. The trial run of 65,000 samples was performed five times for a total of 325,000 samples. The probability of each angle from 0 to 359 is 0.27778 for the exact uniform distribution. Figure 3-13 shows the result of the uniform probability function with two sets of the MC simulation results overlaid onto it. These two scaled sampling arrays show good agreement with the exact uniform distribution; therefore, it can be concluded that the code and the RNG system accurately create the sampling array for the launch angle and functions properly without any logical error. PAGE 77 64 Beta(Degree)Probability (% ) Figure 3-13. Verification of sampling algorithm for rotational angle beta. Two sets of MC trials of 65,000 samples with the uniform probability function. 3.2.4 Verification of the Pathlength Distribution The verification of the pathlength distribution was performed by sampling the launch angle 65,000 times. The trial run of 65,000 samples was performed five times for a total of 325,000 samples. In the process of verification, the mean pathlength was 15.0 m and the standard deviation was 4.0 m based on the assumption in the section 3.1.4. Each value of the mean and standard deviation for the five trials is put on a list in Table 3-2 including the averages of them. The average mean of 14.96 m and standard deviation of 3.89 m correspond with the input vales of the mean and the standard deviation over the pathlength distribution. Table 3-2. The five trials of the mean and standard deviation for sampling the pathlength, along with the average of them. Actual Trial 1 Trial 2 Trial 3 Trial 4 Trial 5 Average Mean ( x ) 15.0 14.98 14.95 14.97 14.96 14.94 14.96 Standard Deviation x 4.0 3.89 3.91 3.89 3.90 3.87 3.89 PAGE 78 65 Figure 3-14 shows the result of the lognormal probability function with two sets of the MC simulation results overlaid onto it. These two scaled sampling arrays show good agreement with the exact lognormal distribution; therefore, it can be concluded that the code and the RNG system accurately create the sampling array for the pathlength and operates properly without any logical error. 0246810126111621263PathlengthProbability(%) 1 Zold Distribution Trial 1 Trial 2 Figure 3-14. Verification of sampling algorithm for pathlength. Two sets of MC trials of 65,000 samples each with the lognormal probability function. 3.2.5 Verification of Photon Progression For verifying of photon progression, the position of individual photons was plotted and tabulated as they traveled through the medium. The tracking code was verified based on the fact that the tracking of a photon terminates when the photon is out of bounds or absorbed. In Table 3-3, an example of the simulation output of the scattering angles, the rotational angle, the pathlength and the absolute position (x, y and z) of a representative photon is tabulated as it progresses through the photon tracking code. In this example, the initial launch angle was 1.13139 with the initial rotational angle of 78.0. The calculated PAGE 79 66 radius was recorded to recognize whether the photon was in the medium. The explanation of how to calculate the radius was introduced in Appendix A. The radius of the final position of the photon was 7482.82 that was less than the inner radius of the cornea (7500m). This indicates that the photon exited through the back surface of the cornea, and the tracking of it ended because the photon was placed out of the boundary. Table 3-3. The absolute coordinates of a photon obtained from MC simulation and independent calculation as it progresses through a simulated human cornea. Monte Carlo Program Output Calculated Values L X Y Z Radius X Y Z 0.00 0.00 0.00 8000.00 29.88 0.12 0.58 29.87 7970.13 0.12 0.58 29.87 2 3 22.72 0.34 1.80 52.56 7947.44 0.34 1.80 52.56 41 115 22.24 -13.52 -1.31 69.67 7930.34 -13.52 -1.31 69.67 2 148 33.22 -33.33 -6.31 95.87 7904.20 -33.33 -6.30 95.87 2 -71 43.72 -59.21 -14.38 130.17 7870.07 -59.21 -14.38 130.17 1 163 40.38 -83.04 -21.17 162.06 7838.41 -83.04 -21.16 162.05 2 -23 31.31 -101.04 -25.59 187.29 7813.41 -101.04 -25.58 187.29 1 3 29.88 -118.01 -29.39 211.59 7789.36 -118.01 -29.37 211.59 1 76 26.06 -133.08 -32.35 232.64 7768.57 -133.08 -32.34 232.64 0 0 36.09 -153.94 -36.46 261.80 7739.82 -153.94 -36.45 261.80 0 0 24.63 -168.18 -39.27 281.70 7720.23 -168.18 -39.25 281.70 0 0 35.13 -188.49 -43.27 310.08 7692.35 -188.49 -43.25 310.09 2 -136 25.11 -202.28 -46.21 330.85 7671.96 -202.28 -46.19 330.86 0 0 26.06 -216.59 -49.26 352.42 7650.81 -216.60 -49.23 352.42 1 148 27.97 -232.30 -52.24 375.37 7628.35 -232.30 -52.21 375.38 0 0 38.47 -253.90 -56.33 406.94 7597.51 -253.90 -56.31 406.95 0 0 36.09 -274.16 -60.18 436.55 7568.66 -274.16 -60.15 436.56 1 -61 26.06 -288.77 -62.51 458.00 7547.79 -288.77 -62.48 458.02 1 -44 29.40 -304.92 -64.82 482.46 7524.00 -304.92 -64.79 482.47 2 150 28.92 -321.63 -67.23 505.95 7501.25 -321.63 -67.20 505.96 2 9 24.15 -336.22 -69.49 525.06 7482.82 -336.22 -69.46 525.07 PAGE 80 67 In addition, the absolute coordinates independently calculated with coordinate transformation are listed to compare the coordinates from the MC simulation with them. MC program output is in excellent agreement (within round-off error) with the independently calculated values; hence, it could be concluded that the photon tracking code was programmed correctly and the RNG worked properly. The progressions of back-and-front exiting photons are respectively described in Figure 3-15 and 3-16. In addition, examples of photons that are absorbed in the medium are depicted in Figure 3-17. 0100200300400500-150-100-50050100150-150-100-50050100150Z-axisX-axis Y-axis Figure 3-15. Tracks of four back-exiting photons. PAGE 81 68 0100200300400500-1500-1000-500050010001500-1000100200Z-axisX-axisY-axis Figure 3-16. Tracks of four front-exiting photons. 0100200300400500-150-100-50050100150-150-100-50050100150Z-axisX-axisY-axis Figure 3-17. Simulation of photons absorbed. PAGE 82 69 3.3 Monte Carlo Simulation The Monte Carlo computer program was utilized to simulate the corneal scattering problem and the polystyrene suspension problem defined with input bounds, pathlength and wavelength. The main purpose of the simulation is to obtain the distribution of where the photons exit the front surface, and then calculate two color ratios at each collection angle corresponding to the probe used in the experiment. 3.3.1 Corneal Simulations In the corneal simulation, the input values were a refractive index of 1.08-0.01i, a particle size of 6.30 m, an 8000 m outer radius, and a 7500 m inner radius. Three types of wavelengths (532 nm, 594 nm, and 670 nm) and six different mean scattering pathlengths (5 m, 10 m, 15 m, 20 m, 25 m, and 30 m) were used in this simulation; thereby, simulations were performed totally 18 times. The purpose of these simulations was to observe the distribution of back scattered photons at the surface of the anterior cornea. In the program, an angle at which each photon exited was recorded every time it emerged from the anterior corneal surface. It took approximately 4 ~ 24 hours for each simulation to be run with 100,000,000 photons that were sufficient to describe a consistent pattern of distribution. Run time was strongly dependent on pathlength used in these simulations. To make sure that the number of photons was enough to accurately predict the exit distribution, the program was run twice with the same parameters, including a 594 nm wavelength and a mean scattering pathlength of 15.0 m. Outputs of the program run are shown in Figure 3-18. The two sets of program output are indiscernible; hence, 100,000,000 samples are enough to elucidate the angular output. PAGE 83 70 05000100001500020000250003000035000400000102030405Exit Alngle(Degree)Relative Frequency 0 trial 1 trial 2 Figure 3-18. Two sets of program output for 670 nm wavelength and 15.0 m mean scattering pathlength based on 100,000,000 trials each. Each launched photon is categorized into four cases. In other words, it would exit through the front, back, sides of the cornea, or be absorbed in the cornea. In a total of 100,000,000 trials, the number of photons corresponding to each of the four cases for three wavelengths and mean scattering pathlengths is listed in Table 3-4. Most photons exited through the back surface of the cornea, with a small percentage of photons exiting through the front surface of the cornea. The mean scattering pathlength implies how dense particles are in the medium. That is, the smaller the pathlength is, the higher the number of scattering particles is in the medium, which results in the more back scattering of the photons as shown in Table 3-4. Among the photons exiting through the front surface of the cornea at the wavelength of 670 nm, 3.16% was detected at 5.0 m mean scattering pathlength whereas only 0.2% was observed at a 30 m mean scattering pathlength. Few photons exited through the sides of the cornea, and approximately 4.2~5.2% of the total photons were absorbed in the cornea. PAGE 84 71 Table 3-4. The relative frequency of 100,000,000 trial photons corresponding to each four case for three wavelength and mean scattering pathlengths. Pathlength 5 um 10 um 532 nm 594 nm 670 nm 532 nm 594 nm 670 nm BACK 93169313 91975436 91589738 94871481 94588757 94535324 FRONT 1943336 2869121 3162489 664935 860327 883164 SIDE 3736 3475 2932 4338 4307 3556 ABSORBED 4883615 5151968 5244841 4459246 4546609 4577956 Pathlength 15 um 20 um 532 nm 594 nm 670 nm 532 nm 594 nm 670 nm BACK 95266974 95128283 95116621 95450669 95360709 95357075 FRONT 394235 480603 481997 279303 331711 328641 SIDE 4251 4139 3847 3754 4183 3685 ABSORBED 4334540 4386975 4397535 4266274 4303397 4310599 Pathlength 25 um 30 um 532 nm 594 nm 670 nm 532 nm 594 nm 670 nm BACK 95559960 95498401 95494552 95640767 95589451 95591400 FRONT 218059 254712 248545 177019 205999 201350 SIDE 3431 3545 3362 3014 3277 3234 ABSORBED 4218550 4243342 4253541 4179200 4201273 4204016 As wavelength decreases, the dimensionless size parameter increases, which is associated with an increase in overall scattering as shown in Table 3-4. As a result, the distributions of back-exiting photons for each wavelength (532, 594, 670 nm) corresponding to each mean scattering pathlength are shown in Figures 3-19 through 3-24. In conclusion, the back scattering response increases as mean scattering pathlengths decrease. Also, it is observed that the difference between the plots of the responses of the two wavelengths increases as the mean scattering pathlength decreases. As expected, this trend agrees with the results of the experiment of PS suspension. As discussed above, the dimensionless size parameter decreases with increasing wavelength, which is associated with a decrease in forward scattering and an increase in backward scattering, or an increase of the number of photon detected at the front surface. PAGE 85 72 0200004000060000800001000001200001400000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-19. Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 5.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength. 0100002000030000400005000060000700000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-20. Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 10.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength. PAGE 86 73 09000180002700036000450000481216202428323640Exit AlngleRelative Frequency 532nm 594nm 670nm Figure 3-21. Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 15.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength. 050001000015000200002500030000350000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-22. Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 20.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength. PAGE 87 74 05000100001500020000250000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-23. Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 25.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength. 050001000015000200000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-24. Distribution of front-exiting photons for each wavelength (532, 594, 670 nm) at a 30.0 m mean scattering pathlength based on 100,000,000 trials at each wavelength. PAGE 88 75 To examine the two-color response further, three exit regions at each 4, 8, and 12 degree corresponding to the solid angle collection area of three collecting fiber optics were selected. As mentioned in the Chapter 2, the diameter of probe area of the 400 m optics is the only 200 m corresponding to an exit angle difference of approximately 1.4 degree. Therefore, exit responses at three selected angles (4, 8, and 12 degree) with 0.5 of those angles were accordingly summed and then the resulting two-color ratio was calculated with the sum of exit responses. These values and two-color ratios are tabulated in Table 3-5. The two-color ratios for each mean scattering pathlength are shown in Figure 3-25 through 3-27. Table 3-5. The relative frequency corresponding to collection areas and two color ratios for mean scattering pathlengths. Pathlength Wavelength Two Color Ratio 532 nm 594 nm 670 nm 532/594 594/670 532/670 5 m 158356 241995 262647 0.65 0.92 0.60 10 m 49620 65063 66345 0.76 0.98 0.75 Optic2 15 m 28521 35527 35359 0.80 1.00 0.81 20 m 19600 24394 23206 0.80 1.05 0.84 25 m 14647 18152 17002 0.81 1.07 0.86 30 m 11911 14187 13886 0.84 1.02 0.86 5 m 180419 319585 374424 0.56 0.85 0.48 10 m 44427 65578 70158 0.68 0.93 0.63 Optic3 15 m 22902 31425 32384 0.73 0.97 0.71 20 m 15344 19978 20352 0.77 0.98 0.75 25 m 11136 14517 14697 0.77 0.99 0.76 30 m 8684 10990 11268 0.79 0.98 0.77 5 m 134094 211592 237169 0.63 0.89 0.57 10 m 37171 55101 60951 0.67 0.90 0.61 Optic4 15 m 19976 26988 28337 0.74 0.95 0.70 20 m 13283 17218 17876 0.77 0.96 0.74 25 m 9706 12415 12204 0.78 1.02 0.80 30 m 7825 9721 9481 0.80 1.03 0.83 PAGE 89 76 . 0.500.600.700.800.901.001.1005101520253035Pathlength(um)Two-Color Ratio 532/594 594/670 532/670 Figure 3-25. The two color ratio for each mean scattering pathlength at fiber optic 2. . 0.400.500.600.700.800.901.001.1005101520253035Pathlength(um)Two-Color Ratio 532/594 594/670 532/670 Figure 3-26. The two color ratio for each mean scattering pathlength at fiber optic 3. PAGE 90 77 0.500.600.700.800.901.001.1005101520253035Pathlength(um)Two-Color Ratio 532/594 594/670 532/670 Figure 3-27. The two color ratio for each mean scattering pathlength at fiber optic 4. The two-color response of the MC simulations were somewhat different from the PS suspension experimental results. As unexpected, two color ratio increased with increasing pathlengths ranging between 5 m to 30 m. The reason is that the pathlength used in the corneal simulation is much different and does not agree with the pathlengths used in the PS suspension experiments. At the more clinically relevant mean scattering pathlengths of 25 and 30 m, it is shown that some of the two-color ratios increased as the scattering concentration increased (i.e. 594 nm to 670 nm ratio at optic 2). It could be suspected that two color ratios would increase, then peak around a mean scattering pathlength of 25 m, and then the ratio would decrease as mean scattering pathlength increases. In other words, a plot of the two color ratio as function of mean scattering pathlength would be parabolic. Clinically, the range of the mean scattering pathlength used in the simulation is difficult to assess for human cornea and may in fact be too small (i.e. too dense). It could be that the two-color ratio increases with increasing scatter PAGE 91 78 concentration in the clinically relevant range of the mean scattering pathlength. To investigate the two-color scattering response further, the polystyrene suspension scattering program was simulated with the Monte Carlo program to provide a direct comparison of the model and experiments. 3.3.2 Polystyrene Suspension Simulations In the PS suspension simulation, the input values were a refractive index of 1.20-0.00i, a particle size of 7.20 m, a 20 mm outer radius of the medium. The simulation was run with 10,000,000 photons that were still enough to describe the distribution of the front-exiting photon correctly. Three absorption coefficients of the PS suspension corresponding to each wavelength were used because these are dependent on wavelengths (4.4e-8 m-1 of 532 nm, 1.6e-7 m-1 of 594 nm, 4.0e-7 m-1 of 670 nm). Note these coefficients are based on water, the scattering medium. First of all, the boundary of the medium should be determined for the simulation. To determine the inner radius used in PS suspension simulation, the preliminary depth test was performed for the inner radius values ranged from 16 mm to 9 mm on the condition that the outer radius was fixed in 20 mm. Ideally, the inner radius value would approach zero; however, it was limited because of computer run-time consideration. The simulation was run with the wavelength of 594 nm and 120.3 m mean scattering pathlength of the C4 suspension. As a result of the preliminary test, the number of back scattering photons converged in approximately 2,300,000, or 23 percentages of total photons at the inner radius 10 mm. The number of photon for each case at the variety of inner radii is tabulated in Table 3-6 and is shown in Figure 3-28. PAGE 92 79 Table 3-6. The relative frequency of 10,000,000 photons for four cases at each inner radius. Inner radius 16 mm 14 mm 13 mm 12 mm 11 mm 10 mm 9 mm At Back 8760651 7654288 6945718 6155737 5295035 4399909 3515646 At Front 1073065 1602902 1822660 2001100 2148558 2258588 2339644 At Side 158414 729767 1216179 1825238 2536392 3319753 4121465 Absorbed 7870 13043 15443 17925 20015 21750 23245 02 *1064 *1066 *1068 *1061 *10781012141618The Number of Detected PhotonInner Radius (mm) 2.3*106 Figure 3-28. The relative frequency of photon exiting through front surface of 10,000,000 trials at each inner radius. According to the results, the number of photon exiting through side surface and being absorbed in the medium increases while the number of photons exiting through back surface decreases. As the inner radius decreases, the thickness of the medium increases, which causes the enhancement of the absorption. On the other hand, the length of the arc of the back boundary decreases with decreasing inner radius, which results in reduction of the probability of back-exiting photon while increasing the probability of PAGE 93 80 side-exiting photon. Probability that a photon traveled near the back surface returns to the front surface would be very low after a given penetration depth. Therefore, it can be concluded that the loss of photon attributable to the assumption that the inner radius is 10 mm would be negligible. With these conditions, the polystyrene suspension scattering environment was simulated for photon tracking. Wavelengths of 532 nm, 594 nm and 670 nm were used at mean scattering pathlengths corresponding to each suspension (27.2 m, 34.4 m, 59 m and 120.3 m respectively). Figure 3-29 through 3-32 show the front-exiting angle distribution of four suspensions (C4, C3, C2 and C1) for three wavelengths. 050000100000150000200000250000048121620242832364Exit AngleRelative Frequency 0 532nm 594nm 670nm Figure 3-29. The front-exiting angle distribution of C4 suspension for three wavelengths. Mean scattering pathlength corresponding to C4 suspension is 120.3 m. PAGE 94 81 0500001000001500002000002500003000003500000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-30. The front-exiting angle distribution of C3 suspension for three wavelengths. Mean scattering pathlength corresponding to C3 suspension is 59 m. 01000002000003000004000005000000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-31. The front-exiting angle distribution of C2 suspension for three wavelengths. Mean scattering pathlength corresponding to C2 suspension is 34.4 m. PAGE 95 82 01000002000003000004000005000006000000481216202428323640Exit AngleRelative Frequency 532nm 594nm 670nm Figure 3-32. The front-exiting angle distribution of C1 suspension for three wavelengths. Mean scattering pathlength corresponding to C1 suspension is 27.2 m. For the C4 suspension, 44.01~59.51% of the launched photons exited through the back boundary, 13.76~22.58% exited through the front boundary (i.e. incident surface), 26.66~33.2% exited through a side boundary, and 0.06~0.51% of photons were absorbed in the medium. For the C1 suspension, 11.2~16.45% of the launched photons exited through the back boundary, 66.4~52.53% exited through the front boundary (i.e. incident surface), 22.09~30.91% exited through a side boundary, and 0.11~0.7% of photons were absorbed in the medium. For suspension C2 and C3, the similar trend was shown. The scattering responses of the PS suspension simulation is listed in Table 3-7. Overall, as the concentration of suspensions increased, the amount of launched photons exiting through the back boundary decreased and the amount of launched photons exiting through the PAGE 96 83 Table 3-7. The scattering responses of the PS suspension simulation. PS Scattering Response Wavelength Suspension Pathlength Fiber optic 2 Fiber optic 3 Fiber optic 4 C1 27.2 269053 217658 206025 532 nm C2 34.4 228212 182316 165379 C3 59.0 158876 118378 99489 C4 120.3 92869 65175 50407 C1 27.2 385530 330777 308177 594 nm C2 34.4 329319 273242 251245 C3 59.0 232562 179733 156454 C4 120.3 142056 106028 83675 C1 27.2 397239 296188 255173 670 nm C2 34.4 348952 251745 217176 C3 59.0 257371 180921 145515 C4 120.3 165974 110655 81312 front boundary increased. The response was similar to that seen in PS suspension experiment and the corneal simulations. As done in section 3.3.1, the two-color ratio was calculated in the collection angles of about 4, 8, and 12 corresponding to an angular collection area of the actual fiber optic probe. The two-color ratio obtained from summing the response over 0.5 to the left and right of the central perpendicular axis was determined to compare to the laboratory experimental results based on the physical position and size of the fiber optic probe. The results of two-color ratios are listed in Table 3-8. The results of the two color ratio at each collection angle are shown in Figure 3-33 through 3-35. PAGE 97 84 Table 3-8. The two-color ratios of the PS suspension simulation. PS Fiber optic 2 Fiber optic 3 Fiber optic 4 Wavelength Suspension Pathlength Ratio Ratio Ratio C1 27.2 0.70 0.66 0.67 532/594 C2 34.4 0.69 0.67 0.66 C3 59.0 0.68 0.66 0.64 C4 120.3 0.65 0.61 0.60 C1 27.2 0.97 1.12 1.21 594/670 C2 34.4 0.94 1.09 1.16 C3 59.0 0.90 0.99 1.08 C4 120.3 0.86 0.96 1.03 C1 27.2 0.68 0.73 0.81 532/670 C2 34.4 0.65 0.72 0.76 C3 59.0 0.62 0.65 0.68 C4 120.3 0.56 0.59 0.62 0.500.600.700.800.901.000.020.040.060.080.0100.0120.0140.0PathlengthTwo Color Rati o 532/594 594/670 532/670 Figure 3-33. The two color ratios at fiber optic 2. PAGE 98 85 0.500.600.700.800.901.001.101.200.020.040.060.080.0100.0120.0140.0PathlengthTwo Color Rati o 532/594 594/670 532/670 Figure 3-34. The two color ratios at fiber optic 3. 0.500.600.700.800.901.001.101.201.300.020.040.060.080.0100.0120.0140.0PathlengthTwo Color Ratio 532/594 594/670 532/670 Figure 3-35. The two color ratios at fiber optic 4. PAGE 99 86 The two color ratio increases with increased scattering concentration (increased light scatter), with the single exception of the 532 nm/594 nm ratio at fiber optic 3 as shown in experimental results. PAGE 100 CHAPTER 4 INTEGRATED RESULTS AND DISCUSSION 4.1 Integrated Results In the PS suspension experiments and simulation, it was observed that at collection angles (4, 8, and 12) the scattering response increased with increasing concentration. Also, the two-color scattering ratios at those angles increased as scattering concentration increased, with the single exception of the 532 nm/594 nm ratio. It is noted that the experimental two-color ratios still include the PE reference values measured to eliminate the variation due to laser power change. To compare experimental absolute results to MC simulation results, the experimental values were corrected with the correction factor as shown in the equation 4-1. '1'2expexpCCRRcorrected (4-1) where, correctedRexp Corrected experimental ratio expR Experimental ratio '1'2CC Correction factor The detail explanation how to obtain the correction factor is in Appendix C. The correction factor is listed in Table 4-1 along with the standard deviation. Table 4-1. The correction factor. Optic 2 Optic 3 Optic 4 594/532 0.700.07 0.670.05 0.650.10 670/594 0.620.04 0.620.06 0.620.06 670/532 0.430.03 0.420.02 0.410.04 87 PAGE 101 88 Also, the corrected experimental results are tabulated in Table 4-2 through 4-4 with uncorrected experimental values and simulation values. Table 4-2. The corrected experimental results (Fiber optic 2). Fiber Optic2 Suspension Pathlength Experiment Corrected Value Simulation Delta(%) C4 120.3 0.810.13 1.160.18 0.65 43.83 C3 59 0.820.05 1.170.07 0.68 41.95 C2 34.4 0.790.05 1.130.07 0.69 38.86 532/ 594 C1 27.2 0.810.05 1.160.07 0.70 39.51 C4 120.3 0.520.03 0.840.05 0.86 2.48 C3 59 0.550.04 0.890.07 0.90 1.43 C2 34.4 0.600.07 0.970.12 0.94 2.87 594/ 670 C1 27.2 0.620.05 0.920.08 0.97 3.00 C4 120.3 0.420.09 0.980.21 0.56 42.67 C3 59 0.450.05 1.050.11 0.62 40.76 C2 34.4 0.470.07 1.090.17 0.65 40.53 532/ 670 C1 27.2 0.500.06 1.160.14 0.68 41.52 Table 4-3. The corrected experimental results (Fiber optic 3). Fiber Optic3 Suspension Pathlength Experiment Corrected Value Simulation Delta(%) C4 120.3 0.790.06 1.180.09 0.61 48.27 C3 59 0.790.03 1.180.05 0.66 44.03 C2 34.4 0.770.03 1.150.04 0.67 41.70 532/ 594 C1 27.2 0.790.03 1.180.03 0.66 44.03 C4 120.3 0.510.04 0.820.07 0.96 14.31 C3 59 0.540.04 0.870.07 0.99 12.02 C2 34.4 0.620.05 1.000.08 1.09 8.26 594/ 670 C1 27.2 0.650.04 1.050.06 1.12 6.39 C4 120.3 0.400.05 0.950.12 0.59 38.05 C3 59 0.420.03 1.000.07 0.65 35.00 C2 34.4 0.470.04 1.120.09 0.72 35.66 532/ 670 C1 27.2 0.510.04 1.210.09 0.73 39.88 PAGE 102 89 Table 4-4. The corrected experimental results (Fiber optic 4). Fiber Optic4 Suspension Pathlength Experiment Corrected Value Simulation Delta(%) C4 120.3 0.770.07 1.180.10 0.60 49.35 C3 59 0.810.03 1.250.04 0.64 48.64 C2 34.4 0.790.02 1.220.03 0.66 45.70 532/ 594 C1 27.2 0.830.09 1.280.14 0.67 47.53 C4 120.3 0.520.05 0.840.08 1.03 18.57 C3 59 0.550.05 0.890.08 1.08 17.86 C2 34.4 0.620.03 1.000.05 1.16 13.79 594/ 670 C1 27.2 0.660.09 1.060.04 1.21 12.02 C4 120.3 0.400.07 0.980.17 0.62 36.45 C3 59 0.450.04 1.100.10 0.68 38.04 C2 34.4 0.490.03 1.200.07 0.76 36.41 532/ 670 C1 27.2 0.550.05 1.340.12 0.81 39.62 Figure 4-1 through 4-9 compares the MC simulations with the PS suspension experiments. Figure 4-1. The two-color ratio values (532 nm/594 nm) at fiber optic 2 as a function of concentration for the MC simulation and the PS suspension experiment. PAGE 103 90 Figure 4-2. The two-color ratio values (594 nm/670 nm) at fiber optic 2 as a function of concentration for the MC simulation and the PS suspension experiment. Figure 4-3. The two-color ratio values (532 nm/670 nm) at fiber optic 2 as a function of concentration for the MC simulation and the PS suspension experiment. PAGE 104 91 Figure 4-4. The two-color ratio values (532 nm/594 nm) at fiber optic 3 as a function of concentration for the MC simulation and the PS suspension experiment. Figure 4-5. The two-color ratio values (594 nm/670 nm) at fiber optic 3 as a function of concentration for the MC simulation and the PS suspension experiment. PAGE 105 92 Figure 4-6. The two-color ratio values (532 nm/670 nm) at fiber optic 3 as a function of concentration for the MC simulation and the PS suspension experiment. Figure 4-7. The two-color ratio values (532 nm/594 nm) at fiber optic 4 as a function of concentration for the MC simulation and the PS suspension experiment. PAGE 106 93 Figure 4-8. The two-color ratio values (594 nm/670 nm) at fiber optic 4 as a function of concentration for the MC simulation and the PS suspension experiment. Figure 4-9. The two-color ratio values (532 nm/670 nm) at fiber optic 4 as a function of concentration for the MC simulation and the PS suspension experiment. PAGE 107 94 The two-color ratios of the MC simulations are overall in excellent agreement with those of the experimental values with regard to the concentration response (slope). Moreover, the absolute values of 594 nm/670 nm two-color ratios even show the excellent agreement at the optic 2. For the 532 nm/594 nm and 532 nm/670 nm ratios, the absolute values of the experiment are around two times larger than those of the simulation. The PS spheres were modeled as perfect dielectrics while wavelength dependent absorption was assigned to the water medium in the simulation. It is suspected that absorption effect on the wavelength of 532 nm would be excessive in the simulation, which results in twice reduction of back scattering. On the assumption of that, the absolute values of the simulation could be in excellent agreement with those of the experimental values. Overall, the results demonstrate a unique color dependence of light scattering under these dense scattering medium conditions because the color ratio changes with the concentration of scattering particles. 4.2 Slope of the Two-Color Ratios The primary goal of this study is to find the most optimal two-color ratio for inventing measuring equipment that can quantify corneal haze with two color ratio technique. To do this, the slope of the two-color ratio with respect to changing concentration was calculated. Figure 4-10 and 4-11 represent the slope values of experiment and simulation of each two color ratio at each fiber optic. The larger the slope is, the smaller the error of response is; therefore, the most optimal two color ratio is 532 nm/670 nm at fiber optic 4, where the slope value is the largest in both the experiment and MC simulation. This is very promising in the context of probe optimization studies. PAGE 108 95 -0.0200.020.040.060.080.10.120.1412345Fiber OpticSlop e 532/594 594/670 532/670 Figure 4-10. The slope values of the two-color ratios as function of fiber optic position corresponding to collection angle in the PS suspension experiment. 00.010.020.030.040.050.060.071234Fiber OpticSlope 5 532/594 594/670 532/670 Figure 4-11. The slope values of the two-color ratios as function of fiber optic position corresponding to collection angle in the PS suspension simulation. PAGE 109 CHAPTER 5 CONCLUSIONS AND RECOMENDATIONS 5.1 Conclusions In the present study, much was learned about the application of the two-color light scattering system in the PS suspension studies and Monte Carlo simulations. Overall, the experimental and modeling studies were consistent in the context of Mie scattering and dense scattering media. To study about more actual three-dimensional problems, several improvements were accomplished for the MC computer model and probe used in the PS experiments. Mainly, the Monte Carlo program was expanded into three dimensions to apply the more complex corneal geometry and to do more accurate corneal modeling. Also, absorption was considered on the scattering particles in the simulated medium. The probability of absorption was around 5% on average, which made the simulation more ideal. Acceptance angles with respect to the photons exiting the front surface were added to simulate the collection angle of the actual fiber optic probe. Acceptance angles ranged between 0 and 40 were divided into the interval of 0.5 degree, so that the total number of acceptance angles is 80, which increased the spatial resolution as compared to previous studies. In terms of run time, the Monte Carlo program was improved by modifying the scattering angle sampling array. An array having only 181 elements was built with cumulative values corresponding to the scattering angle instead of building the probability array having 80,000,000 elements for scattering angle selection, which enabled to save the computer memory and run time. However, the complexity of the 96 PAGE 110 97 program increased as the code was expanded to three-dimensions. Currently, the simulations typically required 4 to 12 hours to process on a 1-GHz PC platform. Based on the results from the simulation, the fiber optic probe was redesigned to collect light in the optimal collection region. Angles of 4, 8, and 12 were selected as collection angles in the fiber optic probe. Then, it was obtained that the two-color scattering ratios as a function of collection angle. Consequently, the most optimal outcome was the ratio of 532 nm/670 nm at the collection angle of 12, which is promising information to design a practical probe that can assess corneal haze. In the PS suspension experiments and simulation, it was observed that at collection angles (4, 8, and 12) the scattering response increased with increasing concentration. Also, the two-color scattering ratios at those angles increased as scattering concentration increased, with the single exception of the 532 nm/594 nm ratio. Overall, the results of the MC simulations are in excellent agreement with those of the experimental values. This demonstrates a unique color dependence of light scattering under these dense scattering medium conditions because the color ratio changes with the concentration of scattering particles. In the corneal simulation, the scattering response was similar to that shown in the PS suspensions simulation. That is, the scattering response increased with increasing the scattering particles concentration (i.e. the pathlength decreased). That same trend continues with these higher values of pathlength. However, the two-color ratio did not increase with increasing corneal haze as expected. It was elucidated with the use of the current multi-wavelength calculations. As the wavelength is getting longer, the dimensionless size parameter decreases, which is associated with an increase in overall PAGE 111 98 scattering. For small collection angles, enhancement of scattering with increasing is consistent with the PS experimental and Monte Carlo simulation results. However, as the scattering concentration and the collection angle are increased, the two-color response becomes a complex combination of enhanced scattering and light penetration to the larger collection angle. This behavior is represented by the curious observation of cross-over of the corneal exit angle distribution curves. This means that at these short pathlengths (high concentrations) the ratio of a shorter wavelength to a longer wavelength will not always increase with increasing scattering particle concentration, as the overall penetration of light becomes important for large angles. This explains why the two-color ratios at the small collection angle increased with decreasing concentration. At the more clinically relevant mean scattering pathlengths of 25 and 30 m, it was shown that some of two color ratios increased as the scattering concentration increased at fiber optic 2 (see Figure 3-28). It could be suspected that the two-color ratio would increase, peak around mean scattering pathlength of 25 m, and then the ratio would decrease as mean scattering pathlength increased; hence, the plot of the two color ratio as function of mean scattering pathlength would be parabolic. Clinically, the range of the mean scattering pathlength used in the simulation is difficult to assess for human cornea and may in fact be too small (i.e. too dense). It could be that the two-color ratio increases with increasing scatter concentration in the clinically relevant range of the mean scattering pathlength. In conclusion, several refinements were made for the computer code and the design of the fiber optic probe in this study. The multi-color scattering system showed promising results for use as a means to quantify and assess corneal tissue haze. The results indicate PAGE 112 99 that the response does not depend on the individual variation by normalizing the absolute signal level. 5.2 Recommendations Based on the future work recommended in the previous study by Coffey, many improvements were made to the computer model and PS experiments in the present study. Additional work is still needed to further develop the multi-color scattering system as a consistent quantitative technique for the haze assessment of the human cornea. The following issues are recommended for future work. As for the fiber optic probe that was used in the experiments, it must be designed to perfectly contact the surface of the scattering medium. The result of failure to do it leads to measure reflective light from the surface, which defects the purpose of the two-color measurement. The work to develop the probe should involve precision alignment. It is recommended that a scattering particle size distribution should be added to the Monte Carlo program to consider the variation of fibroblastic keratocyte size in the future work on this project. The various size of scattering particles in the medium would show more ideal phenomena of the corneal simulation. Also, more complex corneal geometry is needed to apply to simulation with respect to thickness of the cornea that is not constant. In this study, scattering particles in the cornea were assumed that they distributed all over the cornea equally. Ideally, scattering particles concentrate on the ablated area at which fibroblastic keratocytes activate. Probably the most scattering particles concentrate at the anterior cornea surface in PRK while those concentrate at the area dependent on the ablation in LASIK. It is valuable to consider the concentration of scattering particles at a certain area, which may influence the overall scattering phenomena. PAGE 113 100 In terms of program run time, MC computer code should be continuously improved. As the complexity of the program is increased, the need for more concise programming increases. In addition, the MC modeling effects can be improved by better understanding the role of absorption. The PS spheres were modeled as perfect dielectrics, while wavelength dependent absorption was assigned to the water medium. Such effects in the actual experiment and the MC simulations may explain the different degree of agreement for the various two-color ratios. Lastly, the animal experiments can be conducted with new information about the fiber optic probe. In the similar manner done in the previous study, the new study can be performed with 532 nm/670 nm ratios at the most optical collection angle of 12. PAGE 114 APPENDIX A CALCULATION OF RADIUS AND CORNEAL EXIT ANGLE To locate whether a photon is out of bounds, radius was calculated with new absolute coordinates of the photon. Radius is determined in terms of the center of corneal curvature as denoted in Figure A-1. Figure A-1. Geometry of radius and angle . The following equation A-1 explains how to yield radius. 22zranceDistradius , (A-1) where, 22yxanceDist , (A-2) zRzrouter . (A-3) Also, corneal exit angle () at front and back surface was calculated by the equation A-4 for gaining corneal exit angular distributions. )tan(tan1zrceDis . (A-4) 101 PAGE 115 APPENDIX B FORTRAN MONTE CARLO PROGRAM CODE C****************************************************************************** C CALCULATE TRACK OF PARTICLE FOR 3 DIMENSION C****************************************************************************** USE DFPORT USE DFLIB CHARACTER(9) TODAY,NOW,E INTEGER NM(181), ITEST(181),ITEST2(5000) INTEGER ITOTAL,KSTOP,NRUNS,COUNT,FLAG,EXTINCT INTEGER LESS,GREAT,SIDE,A INTEGER*4 ISEED1,K,J,NSUM,NSTOP,NNSUM,NNSTOP,MSTOP,MMSTOP,OSTOP INTEGER*4 NNNSTOP,NNNSUM,MMMSTOP,MMMSUM,NEWTOW,SCALE REAL*8 ITHETA,KTHETA,KBETA,KIBETA REAL*8 PERCENTL,PERCENTG,PERCENTS,PERCENTE REAL*8 THETD(100),ELTRMX(4,100,2),DEG(100), ARRAY(181,3) REAL*8 QEXT,QSCAT,CTBRQS, CUP(181), CNORM(181), CNM(181) REAL*8 W,WAVEL,ANGLE,CTOT,CMIN,CMULT,VCHECK REAL*8 KALPHA,TOW,KABS,CCSUM(181),NSCANG REAL*8 SIG,DMOD,PDMOD,PD,DMAX,DMIN,DELD,XN,ALPHA,DIST(101,2) REAL*8 SIGMA,DMEAN,PAMEAN,PA,AMIN,AMAX,DELA,AN,PATH,DIST2(101,2) REAL*8 ASIG,AMEAN,TEST,TTEST,STEST,PSET(5000),PSET2(5000) REAL*8 BMEAN,BSIG,PBMEAN,BMIN,BMAX,DELB,BN,PB,TTTEST,DDANG REAL*8 DIST3(361,2),PSET3(36100),ESCINDT,TESTAR(181,2) REAL*8 DANG,SCANG,TRAVEL,ESCIND,ESCAR(361,81,4),TESCAR(81,4) REAL*8 RIN,ROUT,YR,RADIUS,DISTANCE,GAMMA,GAMMAH,BETAH,ABETAH REAL*8 R11(2300),R21(2300),R31(2300),R12(2300),R22(2300) REAL*8 R32(2300),R13(2300),R23(2300),R33(2300) OPEN(1,FILE = 'F670-15.OUTPUT') C**** THE START TIME AND DATE OF THE CODE RUN CALL DATE(TODAY) CALL TIME(NOW) WRITE(1,*) TODAY WRITE(1,*) NOW PI = 3.14159 C**** INPUT WAVELENGTH (m) AND REFRACTIVE INDICES: C**** WAVELENGTHS FOR EXPERIMENT (532.0 nm, 594.0 nm,670.0 nm) C WAVEL = 594.0E-9/1.33 C**** CONSOL INPUT ONE OF WAVELENGTHS WRITE(*,*) "INPUT WAVELENGTH(nm) YOU WANT:" READ(*,*) W WAVEL = W*(1.E-9)/1.33 XREAL = 1.44/1.33 XIMAG = 0.01 C**** INPUT DIAMETER (m) DIA = 6300.0E-9 WAVEOUT = WAVEL*1.E9 WRITE(1,100) WAVEOUT WRITE(1,120) XREAL,XIMAG 102 PAGE 116 103 SIZE = DIA*1.E9 WRITE(1,110) SIZE C WRITE(1,130) WRITE(1,101) X = PI*DIA/WAVEL C****************************************************************************** C SET UP PROBABILITY ARRAY (PROBAR) FOR SCATTERING ANGLE THETA C****************************************************************************** C**** SETS VALUES IN THETA AND DEG TO 0 TO 90 JX = 91 ANGLE = 0.0 DO 300 I=1,91 THETD(I) = ANGLE DEG(I) = ANGLE ARRAY(I,1) = ANGLE ANGLE = ANGLE + 1.0 300 CONTINUE C**** CALCULATES CVV AND CHH VALUES FOR ANGLES 0 TO 90 CALL DBMIE(X,XREAL,XIMAG,THETD,JX,QEXT,QSCAT,CTBRQS,ELTRMX) DO 310 I=1,91 CVV1 = (WAVEL*100.)**2*ELTRMX(2,I,1)/(4.*PI**2) CHH1 = (WAVEL*100.)**2*ELTRMX(1,I,1)/(4.*PI**2) ARRAY(I,2) = CVV1 ARRAY(I,3) = CHH1 C WRITE(1,140) DEG(I),CVV1,CHH1 310 CONTINUE C**** SETS VALUES IN DEG TO 91 TO 180 JX = 90 ANGLE = 91.0 DO 301 I=1,90 THETD(I) = 180. ANGLE DEG(I) = ANGLE J=I+91 ARRAY(J,1)=ANGLE ANGLE = ANGLE + 1.0 301 CONTINUE C**** CALCULATES CVV AND CHH VALUES FOR ANGLES 91-180 CALL DBMIE(X,XREAL,XIMAG,THETD,JX,QEXT,QSCAT,CTBRQS,ELTRMX) DO 311 I=1,90 CVV2 = (WAVEL*100.)**2*ELTRMX(2,I,2)/(4.*PI**2) CHH2 = (WAVEL*100.)**2*ELTRMX(1,I,2)/(4.*PI**2) J=I+91 ARRAY(J,2) = CVV2 ARRAY(J,3) = CHH2 C WRITE(1,140) DEG(I),CVV2,CHH2 311 CONTINUE C**** WRITES ARRAY FOR SCATTERING ANGLE THETA DO 320 I=1,181 C WRITE(1,140) ARRAY(I,1), ARRAY(I,2), ARRAY (I,3) 320 CONTINUE C**** CALCULATES Cunpolarized (CUP) AND THE SUM OF CUP (CTOT) CTOT=0 DO 330 I=1,181 CVV3=ARRAY(I,2) CHH3=ARRAY(I,3) CUP(I)=0.5*(CVV3+CHH3) PAGE 117 104 CTOT=CTOT+CUP(I) C WRITE(1,*) CUP(I) 330 CONTINUE C**** NORMALIZES THE CUP ARRAY TO CNORM (Cnormalized) DO 340 I=1,181 CNORM(I)=CUP(I)/CTOT C WRITE(1,*) CNORM(I) 340 CONTINUE C**** CREATES THE SCALED ARRAY CNM (Cnormalized, multiplied) CMIN=MINVAL(CNORM) C WRITE(1,*) CMIN CMULT=1/CMIN CHECK=CMULT*CMIN C WRITE(1,*) CHECK CNM=CMULT*CNORM C**** CONVERTS CNM INTO AN INTEGER ARRAY (NM) AND CUMULATING ARRAY C****(CCSUM) CSUM = 0 DO 350 I=1,181 NM(I)=ANINT(CNM(I)) CSUM=CSUM+NM(I) CCSUM(I)=CSUM C WRITE(1,*) NM(I) 350 CONTINUE C**** BUILD AND INITIALIZE TESTAR ARRAY DO 390 I=1,181 TESTAR(I,1)=I-1 TESTAR(I,2)=0 390 CONTINUE C****************************************************************************** C SET UP UNIFORM PROBABILITY ARRAY (PSET2) FOR LAUNCH ANGLE ALPHA C****************************************************************************** C**** INPUT DATA: AMEAN=MEAN ANGLE ALPHA (DEG) C**** ASIG=STANDARD DEVIATION OF ALPHA (DEG) AMEAN = 0.0 ASIG = 2.33 C**** PARAMETERS: AMIN AND AMAX HAVE 1% PROBABILITY, PAMEAN=PROB AT C****AMEAN PAMEAN=1.0/(ASIG*2.506628) AMIN = AMEAN SQRT(-2.0*(ASIG**2)*LOG(.01)) AMAX = AMEAN + SQRT(-2.0*(ASIG**2)*LOG(.01)) DELA = (AMAX-AMIN)/100.0 AN = 100.0/PAMEAN C**** BUILD DIST2 ARRAY: DIST(I,1)=ALPHA, DIST(I,2)=PROB OF ALPHA ALPHA=AMIN NNSUM=0 DO 400 I=1,101 PA = 1.0/(ASIG*2.506628)*EXP((-(ALPHA-AMEAN)**2)/(2*(ASIG)**2)) DIST2(I,1) = ALPHA DIST2(I,2) = PA*AN NNSUM=NNSUM + ANINT(DIST2(I,2)) C WRITE(1,*) DIST2(I,1), DIST2(I,2) ALPHA=ALPHA + DELA 400 CONTINUE C WRITE(1,*) "NNSUM = ",NNSUM PAGE 118 105 C**** BUILD PROB SET ARRAY: SET SIZE EQUALS TOTAL POPULATION NNSUM K=1 DO 410 I=1,101 NNSTOP=ANINT(DIST2(I,2)) DO 420 J=1,NNSTOP PSET2(K)=DIST2(I,1) K=K+1 420 CONTINUE 410 CONTINUE NNSTOP = K-1 C WRITE(1,*) "PSET2 SIZE =",NNSTOP," LAST SIZE =",PSET2(NNSTOP) C****************************************************************************** C SET UP ZOLD PROBABILITY ARRAY (PSET) FOR PATHLENGTH C****************************************************************************** C**** INPUT DATA 1: SIGMA=STANDARD DEVIATION, Dmean=MEAN DIAMETER (um) C SIGMA = 9.0 C DMEAN = 33.7 C**** CONSOL INPUT DMAN WRITE(*,*) "INPUT MEAN SCATTERING PATHLENGTH(um) YOU WANT:" READ(*,*) DMEAN SIGMA = DMEAN/3.75 C**** PARAMETERS:100 PARTICLES AT MODAL VALUE, Dmin AND Dmax HAVE 1% C**** C**** PROBABILITY C**** SIG=SIGMA NOUGHT, Dmod=MODAL DIAMETER, PDMOD=PROBABILITY AT Dmod SIG = SQRT(LOG((SIGMA/DMEAN)**2+1.0)) DMOD = DMEAN*EXP(-1.50*SIG**2) PDMOD = EXP(-0.50*SIG**2)/(SIG*2.506627*DMOD) DMIN = DMOD*EXP(-1.0*DSQRT(-2.0*SIG**2*LOG(0.01*PDMOD*SIG*DMOD* + 2.506627*EXP(0.5*SIG**2)))) DMAX = DMOD*EXP(1.0*DSQRT(-2.0*SIG**2*LOG(0.01*PDMOD*SIG*DMOD* + 2.506627*EXP(0.5*SIG**2)))) DELD = (DMAX-DMIN)/100.0 XN = 100.0/PDMOD C**** BUILD DIST ARRAY: DIST(I,1)=PATHLENGTH, DIST(I,2)=PROB OF THAT PATHLENGTH PATH = DMIN NSUM = 0 DO 500 I=1,101 PD = EXP(-0.50*SIG**2)*EXP(-0.50*(LOG(PATH/DMOD))**2/SIG**2) + /(SIG*2.506627*DMOD) DIST(I,1) = PATH DIST(I,2) = PD*XN NSUM = NSUM + ANINT(DIST(I,2)) C WRITE(1,*) DIST(I,1),DIST(I,2) PATH = PATH + DELD 500 CONTINUE C WRITE(1,*) "NSUM = ",NSUM C**** BUILD PROB SET ARRAY: SET SIZE EQUALS TOTAL POPULATION, NSUM K = 1 DO 510 I=1,101 NSTOP = ANINT(DIST(I,2)) DO 520 J=1,NSTOP PSET(K) = DIST(I,1) K = K+1 520 CONTINUE 510 CONTINUE PAGE 119 106 NSTOP = K-1 C WRITE(1,*) "PSET SIZE =",NSTOP," LAST SIZE =",PSET(NSTOP) C***************************************************************************** C SET UP NORMAL PROBABILITY ARRAY (PSET3) FOR ANGLE BETA C***************************************************************************** C**** INPUT DATA: BMEAN=MEAN ANGLE BETA (DEG) BMEAN = 0.0 C**** PARAMETERS: BMIN AND BMAX HAVE 1% PROBABILITY, PBMEAN=PROB AT BMEAN PBMEAN=1.0 BMIN = BMEAN 180.0 BMAX = BMEAN + 180.0 DELB = (BMAX-BMIN)/360.0 BN = 100.0/PBMEAN C**** BUILD DIST3 ARRAY: DIST3(I,1)=BETA, DIST3(I,2)=PROBABILITY OF BETA BETA=BMIN NNNSUM=0 DO 600 I=1,361 PB = 1.0 DIST3(I,1) = BETA DIST3(I,2) = PB*BN NNNSUM=NNNSUM + ANINT(DIST3(I,2)) C WRITE(1,*) DIST3(I,1), DIST3(I,2) BETA=BETA + DELB 600 CONTINUE C WRITE(1,*) "NNNSUM = ",NNNSUM C**** BUILD PROB SET ARRAY: SET SIZE EQUALS TOTAL POPULATION NNNSUM K = 1 DO 610 I=1,361 NNNSTOP=ANINT(DIST3(I,2)) DO 620 J=1,NNNSTOP PSET3(K)=DIST3(I,1) K=K+1 620 CONTINUE 610 CONTINUE NNNSTOP = K-1 C WRITE(1,*) "PSET3 SIZE =",NNNSTOP," LAST SIZE =",PSET3(NNNSTOP) C***************************************************************************** C TRACK THE PARTICLES! C***************************************************************************** C**** SET NUMBER OF RUNS(NRUNS) ROUT AND RIN ARE um NRUNS=100000000 ROUT=8000.0 RIN=7500.0 LESS=0 GREAT=0 SIDE=0 EXTINCT=0 C**** RANDOMIZE THE NUMBER GENERATOR: ISEED1 = RND$TIMESEED CALL SEED(ISEED1) C**** BUILD AND INITIALIZE ESCAPE ARRAY C**** FOR BETA FROM 0 TO 360,GAMMA IS FROM 0 TO 40 WITH 0.5 EACH GAPS C DO 697 J=1,361 C WRITE(1,*) 'J=',J PAGE 120 107 C DO 698 I=1,81 C ESCAR(J,I,1)=(I-1)*(50/100.0) C ESCAR(J,I,2)=0.0 C ESCAR(J,I,3)=0.0 C ESCAR(J,I,4)=0.0 C WRITE(1,*) ESCAR(J,I,1), ESCAR(J,I,2), ESCAR(J,I,3), ESCAR(J,I,4) C698 CONTINUE C697 CONTINUE C**** FOR ALL BETA, GAMMA IS ADDED(TESCAR(I,2):TOTAL ESCAPE ARRAY) DO 699 I=1,81 TESCAR(I,1)=(I-1)*(50/100.0) TESCAR(I,2)=0.0 TESCAR(I,3)=0.0 TESCAR(I,4)=0.0 C WRITE(1,*) TESCAR(I,1), TESCAR(I,2), TESCAR(I,3), TESCAR(I,4) 699 CONTINUE C**** BEGIN CALCULATIONS FOR FINDING POSITION OF PARTICLES DO 700 I=1,NRUNS C**** LAUNCH PARTICLE WITH PATHLENGTH, ANGLE ALPHA AND ANGLE BETA C**** NSTOP IS THE SIZE OF PSET(PATHLENGTH) CALL RANDOM(RN) K=ANINT(NSTOP*(1.0-RN)) IF(K.EQ.0) THEN K=1 END IF TRAVEL=PSET(K) C**** NNSTOP IS THE SIZE OF PSET2(LAUNCH ANGLE) CALL RANDOM(RN) K=ANINT(NNSTOP*(1.0-RN)) IF(K.EQ.0) THEN K=1 END IF KALPHA=PSET2(K) IF(KALPHA.LT.0.0) THEN KALPHA=-1*KALPHA END IF C**** NNNSTOP IS THE SIZE OF PSET3(ROTATIONAL ANGLE) IF(KALPHA.EQ.0.0) THEN KIBETA=0.0 ELSE CALL RANDOM(RN) K=ANINT(NNNSTOP*(1.0-RN)) IF(K.EQ.0) THEN K=1 END IF KIBETA=PSET3(K) END IF C WRITE(1,*) "DDANG=",DDANG C**** CALCULATE THE MOVEMENT OF THE PARTICLE X1=TRAVEL*COS(KIBETA*PI/180.0)*SIN(KALPHA*PI/180.0) Y1=TRAVEL*SIN(KIBETA*PI/180.0)*SIN(KALPHA*PI/180.0) Z1=TRAVEL*COS(KALPHA*PI/180.0) C**** FIND THE ABSOLUTE LOCATION OF PARTICLE (ROUT AND RIN ARE um) X=0.0 Y=0.0 Z=0.0 PAGE 121 108 X=X+X1 Y=Y+Y1 Z=Z+Z1 C WRITE(1,*) KALPHA, KIBETA, TRAVEL, X, Y, Z C**** CALCULATE 3*3 ROTATION MATRIX(R) FOR COORDINATE TRANSFORMATION A=1 R11(A)= COS(KIBETA*PI/180.0)*COS(KALPHA*PI/180.0) R12(A)= -SIN(KIBETA*PI/180.0) R13(A)= COS(KIBETA*PI/180.0)*SIN(KALPHA*PI/180.0) R21(A)= SIN(KIBETA*PI/180.0)*COS(KALPHA*PI/180.0) R22(A)= COS(KIBETA*PI/180.0) R23(A)= SIN(KIBETA*PI/180.0)*SIN(KALPHA*PI/180.0) R31(A)= -SIN(KALPHA*PI/180.0) R32(A)= 0.0 R33(A)= COS(KALPHA*PI/180.0) C**** CONTINUE PARTICLE WITH ANGLE THETA,BETA AND PATHLENGTH C**** BEGIN CALCULATIONS, KSTOP IS THE MAXIMUM NUMBER OF ITERATIONS: FLAG=1 COUNT=0 DO WHILE (FLAG.EQ.1) C WRITE(1,*) " IN THE WHILE LOOP" C**** PICK THETA CALL RANDOM(RN) SPOT=INT(CSUM*(1.0-RN)) K=1 IF(SPOT.LE.CCSUM(K).OR.(SPOT.EQ.0.0)) THEN ITHETA=0 ELSE DO 703 J=1,180 IF((SPOT.LE.CCSUM(J+1)).AND.(SPOT.GT.CCSUM(J))) THEN ITHETA=J GO TO 704 END IF 703 CONTINUE END IF 704 KTHETA=ITHETA C WRITE(1,*) KTHETA C**** PICK BETA IF((KTHETA.EQ.0.0).OR.(KTHETA.EQ.180.0)) THEN KBETA=0.0 ELSE CALL RANDOM(RN) K=ANINT(NNNSTOP*(1.0-RN)) IF(K.EQ.0) THEN K=1 END IF KBETA=PSET3(K) END IF C WRITE(1,*) KBETA C**** PICK PATHLENGTH CALL RANDOM(RN) K=ANINT(NSTOP*(1.0-RN)) IF(K.EQ.0) THEN K=1 END IF TRAVEL=PSET(K) PAGE 122 109 C WRITE(1,*) TRAVEL C**** CONSIDERATION OF ABSORPTION EFFECT C**** 95% TRANSMISSION FOR ALL WAVELENGHT FOR HUMEN CORNEA(THICKNESS IS C****500um) C**** KABS AND TRAVEL ARE um KABS=8.595E-5 TOW=EXP(-1*KABS*TRAVEL) C WRITE(1,*) "TOW=",TOW NEWTOW=TOW*1E6 C**** PICK SCALE CALL RANDOM(RN) SCALE=(1E6*(1.0-RN)) C WRITE(1,*) "SCALE=",SCALE IF(SCALE.LT.NEWTOW) THEN GO TO 1000 ELSE IF(SCALE.GE.NEWTOW) THEN C WRITE(1,*) "PARTICLE DISAPPEAR" GO TO 2000 END IF C**** CALCULATE THE MOVEMENT OF THE PARTICLE (RELATVIE COORDINATE SYSTEM) 1000 X2=TRAVEL*COS(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0) Y2=TRAVEL*SIN(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0) Z2=TRAVEL*COS(KTHETA*PI/180.0) C**** FIND THE ABSOLUTE LOCATION OF PARTICLE X=X+R11(A)*X2+R12(A)*Y2+R13(A)*Z2 Y=Y+R21(A)*X2+R22(A)*Y2+R23(A)*Z2 Z=Z+R31(A)*X2+R32(A)*Y2+R33(A)*Z2 C WRITE(1,220) KTHETA, KBETA, TRAVEL, X, Y, Z C**** CALCULATE 3*3 ROTATION MATRIX(R) FOR COORDINATE TRANSFORMATION R11(A+1)=R11(A)*COS(KBETA*PI/180.0)*COS(KTHETA*PI/180.0)+R12(A)* + SIN(KBETA*PI/180.0)*COS(KTHETA*PI/180.0)-R13(A)* + SIN(KTHETA*PI/180.0) R12(A+1)=-R11(A)*SIN(KBETA*PI/180.0)+R12(A)*COS(KBETA*PI/180.0) R13(A+1)=R11(A)*COS(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0)+R12(A)* + SIN(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0)+R13(A)* + COS(KTHETA*PI/180.0) R21(A+1)=R21(A)*COS(KBETA*PI/180.0)*COS(KTHETA*PI/180.0)+R22(A)* + SIN(KBETA*PI/180.0)*COS(KTHETA*PI/180.0)-R23(A)* + SIN(KTHETA*PI/180.0) R22(A+1)=-R21(A)*SIN(KBETA*PI/180.0)+R22(A)*COS(KBETA*PI/180.0) R23(A+1)=R21(A)*COS(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0)+R22(A)* + SIN(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0)+R23(A)* + COS(KTHETA*PI/180.0) R31(A+1)=R31(A)*COS(KBETA*PI/180.0)*COS(KTHETA*PI/180.0)+R32(A)* + SIN(KBETA*PI/180.0)*COS(KTHETA*PI/180.0)-R33(A)* + SIN(KTHETA*PI/180.0) R32(A+1)=-R31(A)*SIN(KBETA*PI/180.0)+R32(A)*COS(KBETA*PI/180.0) R33(A+1)=R31(A)*COS(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0)+R32(A)* + SIN(KBETA*PI/180.0)*SIN(KTHETA*PI/180.0)+R33(A)* + COS(KTHETA*PI/180.0) A=A+1 C**** FIND ANGLE OF THE PARTICLE LOCATION ZR=ROUT-Z DISTANCE=SQRT(X**2+Y**2) C**** FIND THE RADIUS AND ANGLE GAMMA RADIUS=SQRT(DISTANCE**2+ZR**2) PAGE 123 110 GAMMA=ATAN(DISTANCE/ZR)*(180/PI) C**** FOR PARTICLES MOVING OUT OF CORNEA IF(GAMMA.GT.40.03) THEN SIDE=SIDE+1 FLAG=2 C**** FOR PARTICLES MOVING FORWARD(AT CORNEAL BACKSIDE) ELSE IF(RADIUS.LT.RIN) THEN LESS=LESS+1 ESCINDB=INT(80*GAMMA/40)+1 c IF(ESCINDB.EQ.0) THEN c ESCINDB=1 c END IF C**** IN CASE NO CONSIDERATION FOR EACH BETA TESCAR(ESCINDB,2)=TESCAR(ESCINDB,2)+1 c NCOUNT=COUNT+1 c TESCAR(ESCINDB,3)=TESCAR(ESCINDB,3)+NCOUNT C**** IN CASE CONSIDERATION FOR EACH BETA C ESCINDTB=INT(360*BETA/360)+1 C IF(ESCINDTB.EQ.0) THEN C ESCINDTB=1 C END IF C ESCAR(ESCINDTB,ESCINDB,2)=ESCAR(ESCINDTB,ESCINDB,2)+1 FLAG=2 C*********FOR PARTICLES MOVING BACKWARD(AT CORNEAL FRONTSIDE) ELSE IF(RADIUS.GT.ROUT) THEN GREAT=GREAT+1 ESCINDF=INT(80*GAMMA/40)+1 c IF(ESCINDF.EQ.0) THEN c ESCINDF=1 c END IF C**** IN CASE NO CONSIDERATION FOR EACH BETA TESCAR(ESCINDF,3)=TESCAR(ESCINDF,3)+1 NCOUNT=COUNT+1 TESCAR(ESCINDF,4)=TESCAR(ESCINDF,4)+NCOUNT C**** IN CASE CONSIDERATION FOR EACH BETA C ESCINDTF=INT(360*BETA/360)+1 C IF(ESCINDTF.EQ.0) THEN C ESCINDTF=1 C END IF C ESCAR(ESCINDTF,ESCINDF,3)=ESCAR(ESCINDTF,ESCINDF,3)+1 C ESCAR(ESCINDTF,ESCINDF,4)=ESCAR(ESCINDTF,ESCINDF,4)+COUNT FLAG=2 END IF C**** FOR EXTINCT PARTICLES DUE TO ABSORPTION 2000 IF(SCALE.GE.NEWTOW) THEN EXTINCT=EXTINCT+1 FLAG=2 END IF COUNT=COUNT+1 END DO C WRITE(1,*) "COUNT = ", COUNT C WRITE(1,*) X 700 CONTINUE PERCENTL=100*LESS/NRUNS PERCENTG=100*GREAT/NRUNS PRECENTS=100*SIDE/NRUNS PAGE 124 111 PERCENTE=100*EXTINCT/NRUNS C**** THE NUMBER OF PHOTON GONE THROUGH THE CORNEA WRITE(1,*) "BACK = ",LESS," ", "PERCENT=",PERCENTL C**** THE NUMBER OF PHOTON COME BACK TO THE SURFACE OF CORNEA WRITE(1,*) "FRONT = ",GREAT," ", "PERCENT=",PERCENTG C**** THE NUMBER OF PHOTON GONE THROUGH THE CORNEA SIDE WRITE(1,*) "SIDE = ",SIDE," ","PERCENT=",PERCENTS C**** THE NUMBER OF EXTINCT PHOTON DUE TO ABSORPTION WRITE(1,*) "EXTINCT = ",EXTINCT," ","PERCENT=",PERCENTE C**** THE FINAL RESULTS FOR 3D C DO 710 J=1,361 C K=J-1 C WRITE(1,*) "BETA=",K C DO 709 I=1,81 C WRITE(1,221) ESCAR(J,I,1), ESCAR(J,I,2), ESCAR(J,I,3), ESCAR(J,I,4) C709 CONTINUE C710 CONTINUE C**** THE FINAL RESULTS FOR 2D DO 715 I=1,81 WRITE(1,221) TESCAR(I,1),TESCAR(I,2),TESCAR(I,3),TESCAR(I,4) 715 CONTINUE C***************************************************************************** C FORMAT STATEMENTS C***************************************************************************** 100 FORMAT(/1X,'WAVELENGTH (NM) = ',F5.1) 101 FORMAT(/) 110 FORMAT(/1X,'DIAMETER (NM) = ',F6.1) 120 FORMAT(/1X,'REFRACTIVE INDEX = ',F4.2,' ',F4.2,'i') 130 FORMAT(/1X,5x,'ANGLE',8X,'Cvv',12X,'Chh') 140 FORMAT(1x,F10.3,2E15.4) 150 FORMAT(1X,'CMIN = ' E20.11) 170 FORMAT(1X,'CMULT = ' E20.11) 180 FORMAT(1X,'CHECK = ' I20) 190 FORMAT(1X,'CNM = ' E20.11) 200 FORMAT(1X, I10) 210 FORMAT(1X,'PROBAR = ' I10) 220 FORMAT(F7.2,F8.2,F7.2,F7.2,F7.2,F7.2) 221 FORMAT(F7.2,' ',F12.1,F12.1,F12.1) C STOP C**** THE ENDING TIME THE CODE RUN CALL TIME(E) WRITE(1,*) E STOP END C********************************************************************* C SUBROUTINE FOR COMPUTING THE PARAMETERS OF THE ELECTROMAGNETIC C RADIATION SCATTERED BY A SPHERE. THIS SUBROUTINE CARRIES OUT ALL C COMPUTATIONS IN DOUBLE PRECISION ARITHMETIC. C THIS SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE OF C DOWNWARD RECURRENCE RELATIONSHIP. C********************************************************************* SUBROUTINE DBMIE (X,RFR,RFI,THETD,JX,QEXT,QSCAT,CTBRQS,ELTRMX) REAL*8 THETD(100),ELTRMX(4,100,2) REAL*8 QEXT,QSCAT, CTBRQS DIMENSION T(5),TA(4),TB(2),TC(2),TD(2),TE(2) DIMENSION PI(3,91),TAU(3,91),CSTHT(91) PAGE 125 112 DIMENSION SI2THT(91),WFN(2),ACAP(500) COMPLEX RF,RRF,RRFX,WM1,FNA,FNB,TC1,TC2,WFN,ACAP COMPLEX FNAP,FNBP 5 FORMAT(T10,' THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN . 90.0 DEGREES. IT IS ',E15.4) 6 FORMAT (//T10,' PLEASE READ COMMENTS.'//) 7 FORMAT(//T10,'THE VALUE OF THE ARGUMENT JX IS GREATER THAN 100') 8 FORMAT(//T10,'THE UPPER LIMIT FOR ACAPT IS NOT ENOUGH. SUGGEST GET .DETAILED OUTPUT AND MODIFY SUBROUTINE'//) EQUIVALENCE (WFN(1),TA(1)),(FNA,TB(1)),(FNB,TC(1)) EQUIVALENCE(FNAP,TD(1)),(FNBP,TE(1)) IF ( JX .LE. 100 ) GO TO 20 WRITE(6,7) WRITE(6,6) CALL EXIT 20 RF = CMPLX(RFR,-RFI) RRF=1./RF RX=1./X RRFX = RRF * RX T(1) = (X**2)*(RFR**2 + RFI**2) T(1)= SQRT(T(1)) NMX1=1.1*T(1) IF(NMX1.LE.2399) GO TO 21 IX=1 CALL EXIT 21 NMX2 = T(1) IF ( NMX1 .GT. 150 ) GO TO 22 NMX1 = 150 NMX2 = 135 22 ACAP(NMX1+1)=(0.,0.) DO 23 N = 1,NMX1 NN = NMX1 N + 1 ACAP(NN) = (NN+1) * RRFX 1.0 /((NN+1)*RRFX + ACAP(NN+1)) 23 CONTINUE DO 30 J = 1,JX IF ( THETD(J) .LT. 0.0 ) THETD(J) = ABS(THETD(J)) IF ( THETD(J) .GT. 0.0 ) GO TO 24 CSTHT(J) = 1.0 SI2THT(J) = 0.0 GO TO 30 24 IF ( THETD(J) .GE. 90.0 ) GO TO 25 T(1)=(.3141593E1*THETD(J))/180. CSTHT(J) = COS(T(1)) SI2THT(J) = 1.0 CSTHT(J)**2 GO TO 30 25 IF ( THETD(J) .GT. 90.0 ) GO TO 28 CSTHT(J) = 0.0 SI2THT(J) = 1.0 GO TO 30 28 WRITE(6,5) THETD(J) WRITE(6,6) CALL EXIT 30 CONTINUE DO 35 J = 1,JX PI(1,J) = 0.0 PI(2,J) = 1.0 PAGE 126 113 TAU(1,J) = 0.0 TAU(2,J) = CSTHT(J) 35 CONTINUE T(1) = COS(X) T(2) = SIN(X) WM1 = CMPLX( T(1),-T(2)) WFN(1) = CMPLX(T(2),T(1)) WFN(2) = RX * WFN(1) WM1 TC1 = ACAP(1) * RRF + RX TC2 = ACAP(1) * RF + RX FNA = (TC1*TA(3) TA(1))/(TC1*WFN(2) WFN(1)) FNB = ( TC2*TA(3) TA(1))/(TC2 * WFN(2) WFN(1)) N=1 FNAP = FNA FNBP = FNB T(1)=1.5 TB(1) = T(1) * TB(1) TB(2) = T(1) * TB(2) TC(1) = T(1) * TC(1) TC(2) = T(1) * TC(2) DO 60 J = 1,JX ELTRMX(1,J,1) = TB(1) * PI(2,J) + TC(1) * TAU(2,J) ELTRMX(2,J,1) = TB(2) * PI(2,J) + TC(2) * TAU(2,J) ELTRMX(3,J,1) = TC(1) * PI(2,J) + TB(1) * TAU(2,J) ELTRMX(4,J,1) = TC(2) * PI(2,J) + TB(2) * TAU(2,J) ELTRMX(1,J,2) = TB(1) * PI(2,J) TC(1) * TAU(2,J) ELTRMX(2,J,2) = TB(2) * PI(2,J) TC(2) * TAU(2,J) ELTRMX(3,J,2) = TC(1) * PI(2,J) TB(1) * TAU(2,J) ELTRMX(4,J,2) = TC(2) * PI(2,J) TB(2) * TAU(2,J) 60 CONTINUE QEXT=2.*(TB(1)+TC(1)) QSCAT =(TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2 )/.75 CTBRQS = 0. N = 2 65 T(1) = 2*N 1 T(2) = N 1 T(3) = 2 * N + 1 DO 70 J = 1,JX PI(3,J) = (T(1)*PI(2,J)*CSTHT(J)-N*PI(1,J))/T(2) TAU(3,J) = CSTHT(J)*(PI(3,J)-PI(1,J))-T(1)*SI2THT(J)*PI(2,J)+ + TAU(1,J) 70 CONTINUE WM1 = WFN(1) WFN(1) = WFN(2) WFN(2) = T(1)*RX*WFN(1) WM1 TC1 = ACAP(N)*RRF + N*RX TC2 = ACAP(N)*RF + N*RX FNA = (TC1*TA(3)-TA(1))/(TC1*WFN(2) WFN(1)) FNB = (TC2*TA(3)-TA(1))/(TC2*WFN(2) WFN(1)) T(5) = N T(4) = T(1)/(T(5)*T(2)) T(2) = (T(2)*(T(5) + 1.0 ))/T(5) CTBRQS = CTBRQS + T(2)*(TD(1)*TB(1) + TD(2)*TB(2) + TE(1)*TC(1) + + TE(2)*TC(2)) + T(4)*(TD(1)*TE(1) + TD(2)*TE(2)) QEXT = QEXT + T(3)*(TB(1)+TC(1)) T(4) = TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2 PAGE 127 114 QSCAT = QSCAT + T(3) *T(4) T(2) = N*(N+1) T(1) = T(3)/T(2) K = (N/2)*2 DO 80 J = 1,JX ELTRMX(1,J,1)=ELTRMX(1,J,1)+T(1)*(TB(1)*PI(3,J)+TC(1)*TAU(3,J)) ELTRMX(2,J,1)=ELTRMX(2,J,1)+T(1)*(TB(2)*PI(3,J)+TC(2)*TAU(3,J)) ELTRMX(3,J,1)=ELTRMX(3,J,1)+T(1)*(TC(1)*PI(3,J)+TB(1)*TAU(3,J)) ELTRMX(4,J,1)=ELTRMX(4,J,1)+T(1)*(TC(2)*PI(3,J)+TB(2)*TAU(3,J)) IF ( K .EQ. N ) GO TO 75 ELTRMX(1,J,2)=ELTRMX(1,J,2)+T(1)*(TB(1)*PI(3,J)-TC(1)*TAU(3,J)) ELTRMX(2,J,2)=ELTRMX(2,J,2)+T(1)*(TB(2)*PI(3,J)-TC(2)*TAU(3,J)) ELTRMX(3,J,2)=ELTRMX(3,J,2)+T(1)*(TC(1)*PI(3,J)-TB(1)*TAU(3,J)) ELTRMX(4,J,2)=ELTRMX(4,J,2)+T(1)*(TC(2)*PI(3,J)-TB(2)*TAU(3,J)) GO TO 80 75 ELTRMX(1,J,2)=ELTRMX(1,J,2)+T(1)*(-TB(1)*PI(3,J)+TC(1)*TAU(3,J)) ELTRMX(2,J,2)=ELTRMX(2,J,2)+T(1)*(-TB(2)*PI(3,J)+TC(2)*TAU(3,J)) ELTRMX(3,J,2)=ELTRMX(3,J,2)+T(1)*(-TC(1)*PI(3,J)+TB(1)*TAU(3,J)) ELTRMX(4,J,2)=ELTRMX(4,J,2)+T(1)*(-TC(2)*PI(3,J)+TB(2)*TAU(3,J)) 80 CONTINUE IF( T(4) .LT. 1.0E-8 ) GO TO 100 N = N + 1 DO 90 J = 1,JX PI(1,J) = PI(2,J) PI(2,J) = PI(3,J) TAU(1,J) = TAU(2,J) TAU(2,J) = TAU(3,J) 90 CONTINUE FNAP = FNA FNBP = FNB IF ( N .LE. NMX2 ) GO TO 65 WRITE(6,8 ) CALL EXIT 100 DO 120 J = 1, JX DO 120 K = 1,2 DO 115 I = 1,4 T(I) = ELTRMX(I,J,K) 115 CONTINUE ELTRMX(2,J,K) = T(1)**2 + T(2)**2 ELTRMX(1,J,K) = T(3)**2 + T(4)**2 ELTRMX(3,J,K) = T(1)*T(3) + T(2)*T(4) ELTRMX(4,J,K) = T(2)*T(3) T(4)*T(1) 120 CONTINUE T(1)=2.*RX**2. QEXT = QEXT * T(1) QSCAT = QSCAT * T(1) CTBRQS=2.*CTBRQS*T(1) RETURN END C************************************************************** PAGE 128 APPENDIX C PS EXPERIMENTAL DATA CORRECTION In the experiment, the scattering signal is defined as CIS (C-1) Where, S Signal output I Incident energy C Scattering parameter (f(size of scattering particle, number density)) Efficiency of photons entering and exiting the medium Efficiency of fiber/ spectrometer Two color ratios of PS suspensions and PE reference are respectively defined as 2222111121 CICISSRPS , (C-2) '2'2'2'2'1'1'1'121CICISSRPE . (C-3) In order to eliminate the day-to day change of the laser power, the two-color ratios of the PS suspensions were normalized with the ratios of the PE reference values. The two-color ratio obtained in the experimental study is defined as '1'221'1'221'1'221'1'221expCCCCIIIIRRRPEPS (C-4) Based on some assumptions following below, for a single experiment, the laser power is consistent; hence, '11II ;. (C-5) '22II The refractive index is assumed constant over 1 to 2; hence, 115 PAGE 129 116 1'1'221 . (C-6) By definition, '11 ;. (C-7) '22 Hence, the equation C-4 is simplified as '1'221expCCCCR (C-8) In the simulation, 111NCE where, 1N The number of photon launched from the illumination optic. 1E The number of photon detected by the collection optic. 221121NCNCEERMC (C-9) The number of photons launched from the illumination optic is always constant for the simulation. In order to compare Rexp with RMC, it is needed to know '2'1C Cthat can be calculated through the calibration measurement with a power meter. '2'2'2'2'1'1'1'1'2'1CICISSRcal (C-10) In this equation, 1'1'221 ; hence, the equation C-10 is written as 21'2'1'1'2'1'2IISSCC (C-11) PAGE 130 117 '2'1'1'2IISS is calculated with the values measured with the power meter and signal output, and 21 is calculated with calibrated blackbody source, 122121 BBBBIISS (C-12) where and are the signals recorded with the blackbody, and and are the known blackbody emission. 1BS 2BS 1BI 2BI Table C-1. Output of a trial for calculating correction factor. Laser Power Ratio 532 nm 4.1 I532/I594 0.37 594 nm 11.0 I594/I670 1.15 670 nm 9.5 I532/I670 0.43 Signal of PE Optic2 Optic3 Optic4 S594/S532 1.13 1.14 1.14 S670/S594 0.33 0.33 0.37 S670/S532 0.37 0.37 0.42 Signal of BlackBody Optic2 Optic3 Optic4 S594/S532 0.99 1.00 1.04 S670/S594 1.10 1.08 1.08 S670/S532 1.09 1.08 1.12 Incident Energy Optic2 Optic3 Optic4 I532/I594 1.5727 1.4750 1.4570 I594/I670 1.5936 1.7220 1.5140 I532/I670 2.5062 2.5400 2.2060 Efficiency Optic2 Optic3 Optic4 532/594 1.5569 1.4749 1.5160 594/670 1.7530 1.8600 1.6350 532/670 2.7318 2.7434 2.4710 Correction Factor Optic2 Optic3 Optic4 C594/C532 0.65 0.62 0.64 C670/C594 0.67 0.71 0.70 C670/C532 0.43 0.43 0.44 Hence, the correction factor is defined as PAGE 131 118 1221'2'1'1'2'1'2BBBBIISSIISSCC . (C-13) An output of a trial for calculating correction factor is shown in Table C-1. Consequently, new corrected two-color ratio is obtained by dividing two-color ratio from the PS experiments by the calculated correction factor. '1'2expexpCCRRcorrected . (C-14) PAGE 132 LIST OF REFERENCES Al-qahtani J. M., Mclean I. W., Weiblinger R. P. and Ediger M. N., 2001, Preliminary in vitro study of the histological dffects of low fluence 193 nm excimer laser irradiation of corneal tissue. 17, 105-109. Andrade H. A., McDonald M. B., Liu J. C., Abdelmegeed M., Varnell R. and Sunderland G., 1990, Evaluation of an opacity lensometer for determining corneal clarity following excimer laser photoablation. Refractive & Corneal Surgery. 6:346-351. Bevilacqua F., Piguet D., Marquet P., Gross J.D., Tromberg B.J. and Depeursinge C., 1999, In vivo local determination of tissue optical properties: applications to human brain. Applied Optics. 38(22), 4939-4949. Braunstein R. E., Jain S., McCally R. L., Stark W. J., Connolly P. J. and Azar D. T., 1996, Objective measurement of corneal light scattering after excimer laser keratectomy. Ophthalmology. 103(2), 439-443. Coffey L.M., 2001, An investigation of light scattering as a means to monitor corneal haze. Masterâ€™s. Thesis, Department of Mechanical and Aerospace Engineering, University of Florida. Crane III C.D. and Duffy J., 1998, Kinematic analysis of robot manipulators, University of Florida, Chapter 2. Dave J. V., 1968, Subroutines for computing the parameters of the electromagnetic radiation scattered by a sphere. IBM Palo Alto Scientific Center Report No. 320-3237. del Val J. A., Barrero S., Yanez B., Merayo J., Aparicio J. A., Gonzalez V. R., Pastor J. C. and Mar S., 2001, Experimental measurement of corneal haze after excimer laser keratectomy. Applied Optics. 40(10), 1727-1734. Eiferman R. A., Oâ€™Neill K. P., Forgey D. R. and Cook Y. D., 1991, Excimer laser photorefractive keratectomy for myopia: six-month results. Refractive & Corneal Surgery. 7, 344-347. Espenscheid W. F., Kerker M. and Matijevic E., 1964, Logarithmic distribution function for colloidal particles. Journal of Physical Chemistry. 68, 3093-3097. 119 PAGE 133 120 Fantes F. E., Hanna K. D., Waring III G. O., Pouliquen Y., Thompson K. P. and Savoldelli, M., 1990, Wound healing after excimer laser keratomileusis (photorefractive keratectomy) in monkeys. Archives in Ophthalmology. 108, 665-675. Gartry D. S., Muir M. G. K. and Marshall J., 1991, Photorefractive keratectomy with an argon fluoride excimer laser: a clinical study. Refractive & Corneal Surgery. 7, 420-435. Hahn D. W., Flower W. L. and Hencken K. R., 1997, Discrete particle detection and metal emissions monitoring using laser-induced breakdown spectroscopy. Applied Spectroscopy. 51(12), 1836-1844. Jain S., Khoury J. M., Chamon W. and Azar D. T., 1995, Corneal light scattering after laser in situ keratomileusis and photorefractive keratectomy. American Journal of Ophthalmology. 120(4), 532-534. Jester J. V., Moller-Pedersen T., Huang J., Sax C. M., Kays W. T., Cavanagh H. D., Petroll W. M. and Piatigorsky J., 1999a, The cellular basis of corneal transparency: Evidence for â€œcorneal crystallins.â€ Journal of Cell Science. 112, 613-622. Jester J. V., Petroll W. M. and Cavanagh H. D., 1999b, Corneal stromal wound healing in refractive surgery: the role of myofibroblasts. Progress in Retinal and Eye Research. 18(3), 311-356. Jones, A. R., 1979, Scattering of electromagnetic radiation in particulate laden fluids. Prog. Energy Combust. Sci.. 5, 73-96. Kaji Y., Obata H., Usui T., Soya K., Machinami R., Tsuru T. and Yamashita H., 1998, Three-dimensional organization of collagen fibrils during corneal stromal wound healing after excimer laser keratectomy. Journal of Cataract & Refractive Surgery. 24, 1441-1446. Kaplan B., Ledanois G. and Drevillon B., 2001, Mueller matrix of dense polystyrene latex sphere suspensions: Measurements and Monte Carlo simulation. Applied Optics. 40(16), 2769-2777. Kerker M., 1969, The Scattering of Light and Other Electromagnetic Radiation. Academic, New York. Key H., Davies E. R., Jackson P. C. and Wells P. N. T., 1991, Monte Carlo modeling of light propagation in breast tissue. Physics in Medicine and Biology. 36(5), 591-602. Kienle A., Patterson M. S., Ott L. and Steiner R., 1996, Determination of the scattering coefficient and the anisotropy factor from laser Doppler spectra of liquids including blood. Applied Optics. 35(19), 3404-3412. PAGE 134 121 Klyce S. D. and Beuerman R. W., 1998, Structure and function of the cornea. The Cornea: Second Edition (Edited by Kaufman, H. E., Barron, B. A. and McDonald, M. B.), pp. 3-50. Butterworth-Heinemann, Boston. Koelink M. H., de Mul F. F. M., Greve J., Graaff R., Dassel A. C. M. and Aarnoudse J. G., 1994, Laser Doppler blood flowmetry using two wavelengths: Monte Carlo simulations and measurements. Applied Optics. 33(16), 3549-3558. Lohmann C. P., Timberlake G. T., Fitzke F. W., Gartry D. S., Muir M. K. and Marshall J., 1992, Corneal light scattering after excimer laser photorefractive keratectomy: the objective measurements of haze. Refractive & Corneal Surgery. 8, 114-121. Manche E. E., Carr J. D., Haw W. W. and Hersh P. S., 1998, Excimer laser refractive surgery. The Western Journal of Medicine. 169(1), 30-38. Martini F. H., 2001, Anatomy & Physiology. Prentice Hall, Upper Saddle River. McDonald M. B. and Chitkara D., 1998, Principles of excimer laser photoablation. The Cornea: Second Edition (Edited by Kaufman, H. E., Barron, B. A. and McDonald, M. B.), pp. 973-980. Butterworth-Heinemann, Boston. Moller-Pedersen T., Cavanagh H. D., Petroll W. M. and Jester J. V., 2000, Stromal wound healing explains refractive instability and haze development after photorefractive keratectomy. Ophthalmology. 107(7),1235-1245. Moller-Pedersen T., Voel M., Li H. F., Petroll M., Cavanagh D. and Jester J. V., 1997, Quantification of stromal thinning, epithelial thickness, and corneal haze after photorefractive keratectomy using in vivo confocal microscopy. Ophthalmology. 104, 360-368. Niemz M. H., 1996, Laser-Tissue Interactions. Springer-Verlag, New York. Odrich M. G. and Greenberg K. A., 1998, Photorefractive keratectomy for myopia. The Cornea: Second Edition (Edited by Kaufman, H. E., Barron, B. A. and McDonald, M. B.), pp. 981-997. Butterworth-Heinemann, Boston. Olsen T., 1982, Light scattering from the human cornea. Invest Ophathalmol Vis Sci. 23, 81-86 Parrish C. M. and Chandler J. W., 1998, Corneal trauma. The Cornea: Second Edition (Edited by Kaufman, H. E., Barron, B. A. and McDonald, M. B.), pp. 633-671. Butterworth-Heinemann, Boston. Pogue B. W. and Burke G., 1998, Fiber-optic bundle design for quantitative fluorescence measurement from tissue. Applied Optics. 37(31), 7429-7436. PAGE 135 122 Rouweyha R. M., Chuang A. Z., Mitra S., Philips C. B. and Yee R. W. 2002, Laser epithelial keratomileusis for myopia with the autonomous laser. Journal of Refractive Surgery. Volume 18. Sacks Z. S., Craig D. L., Kurtz R. M., Juhasz T. and Mourou G., 1998, Spatially resolved transmission of highly focused beams through cornea and sclera between 1400 and 1800 nm. Saratov Conference, Department of Ophthalmology, University of Michigan., http://www/eecs/umich.edu/CUOS/index.html . July/7/2002 Seiler T., Kahle G. and Kriegerowski M., 1990, Excimer laser (193 nm) myopic keratomileusis in sighted and blind human eyes. Refractive & Corneal Surgery. 6, 165-173. Smith G. T. H., Brown N. A. P. and Shun-Shin A., 1990, Light scatter from the central human cornea. Eye. 4, 584-588. PAGE 136 BIOGRAPHICAL SKETCH Kibum Kim, the first son of Yong-Hyun Kim and Pill-Young Kimâ€™s two sons, was born in Dae-Jeon, Republic of Korea, on February 4, 1975. He grew up and spent most of his early life in the city. In 2000, he graduated and achieved the bachelorâ€™s degree in naval architecture and ocean Engineering from the Chung-Nam National University in the Republic of Korea. In his senior year, he built A human-powered ship with water foils in Korea for the first time. He then obtained the opportunity of having experience about American culture and studying English at the University of Alabama during the following year. In August of 2001, he began pursuing a Master of Science degree in mechanical engineering at the University of Florida, which incorporated the efforts of this masterâ€™s thesis. Kibum will keep going on pursuing his Ph.D. in mechanical engineering at the same university beginning August 2003. 123 |