forked from ankrh/MCmatlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Example1_StandardTissue.m
108 lines (93 loc) · 5.62 KB
/
Example1_StandardTissue.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
%% Description
% In this introductory example, a block of "standard tissue" (mu_a = 1,
% mu_s = 100, g = 0.9) is illuminated by a pencil beam (infinitely thin
% beam). A small slice of air is present in the top of the simulation
% volume. nx and ny are set to odd values so that when the pencil beam is
% launched at x = y = 0 and travels straight down, it travels along a
% well-defined center column of voxels (the middle of the 51st column). Use
% the log10 plot checkbox in the visualizations to better see the fluence
% rate and absorption distribution in the MC result.
%% 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 = 101; % Number of bins in the x direction
model.G.ny = 101; % Number of bins in the y direction
model.G.nz = 150; % 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 = .15; % [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.simulationTimeRequested = .1; % [min] Time duration of the simulation
model.MC.matchedInterfaces = true; % Assumes all refractive indices are the same
model.MC.boundaryType = 1; % 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 = 532; % [nm] Excitation wavelength, used for determination of optical properties for excitation light
model.MC.lightSource.sourceType = 0; % 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)
% For a pencil beam, the "focus" is just a point that the beam goes
% through, here set to be the center of the cuboid:
model.MC.lightSource.xFocus = 0; % [cm] x position of focus
model.MC.lightSource.yFocus = 0; % [cm] y position of focus
model.MC.lightSource.zFocus = model.G.Lz/2; % [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
% These lines will run the Monte Carlo simulation with the provided
% parameters and subsequently plot the results:
model = runMonteCarlo(model);
model = plot(model,'MC');
%% Geometry function(s) (see readme for details)
% A geometry function takes as input X,Y,Z matrices as returned by the
% "ndgrid" MATLAB function as well as any parameters the user may have
% provided in the definition of Ginput. It returns the media matrix M,
% containing numerical values indicating the media type (as defined in
% mediaPropertiesFunc) at each voxel location.
function M = geometryDefinition(X,Y,Z,parameters)
zSurface = 0.03;
M = ones(size(X)); % Air
M(Z > zSurface) = 2; % "Standard" tissue
end
%% Media Properties function (see readme for details)
% The media properties function defines all the optical and thermal
% properties of the media involved by filling out and returning a
% "mediaProperties" array of "mediumProperties" objects with various
% properties. The j indices are the numbers that are referred to in the
% geometry function (in this case, 1 for "air" and 2 for "standard tissue")
% See the readme file or the examples for a list of properties you may
% specify. Most properties may be specified as a numeric constant or as
% function handles.
%
% The function must take one input; the cell array containing any special
% parameters you might specify above in the model file, for example
% parameters that you might loop over in a for loop. In most simulations
% this "parameters" cell array is empty. Dependence on wavelength is shown
% in examples 4 and 23. Dependence on excitation fluence rate FR,
% temperature T or fractional heat damage FD can be specified as in
% examples 12-15.
function mediaProperties = mediaPropertiesFunc(parameters)
% Always leave the following line in place to initialize the
% mediaProperties array:
mediaProperties = MCmatlab.mediumProperties;
% Put in your own media property definitions below:
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