forked from ankrh/MCmatlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathExample8_CollectedFluenceRate.m
119 lines (99 loc) · 6.86 KB
/
Example8_CollectedFluenceRate.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
%% Description
% This example shows how to use one of the "deposition criteria" to
% calculate and plot the fluence rate (and the other outputs) of only that
% light which ends up on the light collector. The geometry is a block of
% standard tissue, into which we launch a Gaussian beam. We have a light
% collector looking at y = -0.02.
%
% We specify that we only want to look at those photons that end up on the
% light collector with model.MC.depositionCriteria.onlyCollected = true;
% See examples 19 and 25 for more use of the deposition criteria.
%
% In the fluence rate plot, you can see how the photons all start at the
% source and end at the light collector. This is also very clear in the
% photon paths plot.
%
% To use a light collector, the cuboid boundary type towards the detector
% has to be set to "escaping". Additionally, the voxels touching that
% boundary must have a refractive index of 1.
%
% This example also shows two other features: (1) That the Monte Carlo
% simulation can be set to launch a set number of photons rather than run
% for a set time using the MC.nPhotonsRequested property, and (2) that by
% setting the boundaryType flags to 2, we can allow photons to travel
% outside the cuboid in the x, y, and +z directions as seen in the photon
% paths illustrations, although absorption and fluence rate data is still
% only tracked within the main cuboid.
%% MCmatlab abbreviations
% G: Geometry, MC: Monte Carlo, FMC: Fluorescence Monte Carlo, HS: Heat
% simulation, M: Media array, FR: Fluence rate, FD: Fractional damage.
%
% There are also some optional abbreviations you can use when referencing
% object/variable names: LS = lightSource, LC = lightCollector, FPID =
% focalPlaneIntensityDistribution, AID = angularIntensityDistribution, NI =
% normalizedIrradiance, NFR = normalizedFluenceRate.
%
% For example, "model.MC.LS.FPID.radialDistr" is the same as
% "model.MC.lightSource.focalPlaneIntensityDistribution.radialDistr"
%% Geometry definition
MCmatlab.closeMCmatlabFigures();
model = MCmatlab.model;
model.G.nx = 100; % Number of bins in the x direction
model.G.ny = 100; % Number of bins in the y direction
model.G.nz = 100; % Number of bins in the z direction
model.G.Lx = .1; % [cm] x size of simulation cuboid
model.G.Ly = .1; % [cm] y size of simulation cuboid
model.G.Lz = .1; % [cm] z size of simulation cuboid
model.G.mediaPropertiesFunc = @mediaPropertiesFunc; % Media properties defined as a function at the end of this file
model.G.geomFunc = @geometryDefinition; % Function to use for defining the distribution of media in the cuboid. Defined at the end of this m file.
model = plot(model,'G');
%% Monte Carlo simulation
model.MC.nPhotonsRequested = 5e6; % # of photons to launch
model.MC.nExamplePaths = 20;
model.MC.matchedInterfaces = true; % Assumes all refractive indices are the same
model.MC.boundaryType = 2; % 0: No escaping boundaries, 1: All cuboid boundaries are escaping, 2: Top cuboid boundary only is escaping, 3: Top and bottom boundaries are escaping, while the side boundaries are cyclic
model.MC.wavelength = 400; % [nm] Excitation wavelength, used for determination of optical properties for excitation light
model.MC.lightSource.sourceType = 4; % 0: Pencil beam, 1: Isotropically emitting line or point source, 2: Infinite plane wave, 3: Laguerre-Gaussian LG01 beam, 4: Radial-factorizable beam (e.g., a Gaussian beam), 5: X/Y factorizable beam (e.g., a rectangular LED emitter)
model.MC.lightSource.focalPlaneIntensityDistribution.radialDistr = 1; % Radial focal plane intensity distribution - 0: Top-hat, 1: Gaussian, Array: Custom. Doesn't need to be normalized.
model.MC.lightSource.focalPlaneIntensityDistribution.radialWidth = .01; % [cm] Radial focal plane 1/e^2 radius if top-hat or Gaussian or half-width of the full distribution if custom
model.MC.lightSource.angularIntensityDistribution.radialDistr = 1; % Radial angular intensity distribution - 0: Top-hat, 1: Gaussian, 2: Cosine (Lambertian), Array: Custom. Doesn't need to be normalized.
model.MC.lightSource.angularIntensityDistribution.radialWidth = 0; % [rad] Radial angular 1/e^2 half-angle if top-hat or Gaussian or half-angle of the full distribution if custom. For a diffraction limited Gaussian beam, this should be set to model.MC.wavelength*1e-9/(pi*model.MC.lightSource.focalPlaneIntensityDistribution.radialWidth*1e-2))
model.MC.lightSource.xFocus = 0; % [cm] x position of focus
model.MC.lightSource.yFocus = 0.025; % [cm] y position of focus
model.MC.lightSource.zFocus = 0; % [cm] z position of focus
model.MC.lightSource.theta = 0; % [rad] Polar angle of beam center axis
model.MC.lightSource.phi = 0; % [rad] Azimuthal angle of beam center axis
model.MC.useLightCollector = true;
model.MC.lightCollector.x = 0; % [cm] x position of either the center of the objective lens focal plane or the fiber tip
model.MC.lightCollector.y = -0.025; % [cm] y position
model.MC.lightCollector.z = 0; % [cm] z position
model.MC.lightCollector.theta = 0; % [rad] Polar angle of direction the light collector is facing
model.MC.lightCollector.phi = pi/2; % [rad] Azimuthal angle of direction the light collector is facing
model.MC.lightCollector.f = .1; % [cm] Focal length of the objective lens (if light collector is a fiber, set this to Inf).
model.MC.lightCollector.diam = .1; % [cm] Diameter of the light collector aperture. For an ideal thin lens, this is 2*f*tan(asin(NA)).
model.MC.lightCollector.fieldSize = .04; % [cm] Field Size of the imaging system (diameter of area in object plane that gets imaged). Only used for finite f.
model.MC.lightCollector.NA = 0.22; % [-] Fiber NA. Only used for infinite f.
model.MC.lightCollector.res = 50; % X and Y resolution of light collector in pixels, only used for finite f
model.MC.depositionCriteria.onlyCollected = true;
model = runMonteCarlo(model);
model = plot(model,'MC');
%% Geometry function(s) (see readme for details)
function M = geometryDefinition(X,Y,Z,parameters)
zSurface = 0.01;
M = ones(size(X)); % Air
M(Z > zSurface) = 2; % "Standard" tissue
end
%% Media Properties function (see readme for details)
function mediaProperties = mediaPropertiesFunc(parameters)
mediaProperties = MCmatlab.mediumProperties;
j=1;
mediaProperties(j).name = 'air';
mediaProperties(j).mua = 1e-8; % Absorption coefficient [cm^-1]
mediaProperties(j).mus = 1e-8; % Scattering coefficient [cm^-1]
mediaProperties(j).g = 1; % Henyey-Greenstein scattering anisotropy
j=2;
mediaProperties(j).name = 'standard tissue';
mediaProperties(j).mua = 1; % Absorption coefficient [cm^-1]
mediaProperties(j).mus = 100; % Scattering coefficient [cm^-1]
mediaProperties(j).g = 0.9; % Henyey-Greenstein scattering anisotropy
end