forked from ankrh/MCmatlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Example7_FarField.m
129 lines (111 loc) · 7.45 KB
/
Example7_FarField.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
120
121
122
123
124
125
126
127
128
129
%% Description
% Here, we demonstrate calculation and plotting of the far field angular
% distribution of the "escaped" photons, both for excitation light and for
% fluorescence light. It is enabled by specifying the optional
% "farFieldRes" parameter, which serves also to specify the resolution you
% want of the far field distribution. The geometry is almost the same
% fluorescing cylinder as example 5, but now illuminated by a tilted
% Gaussian beam and now with only one type of fluorescer.
%
% In the far field of the excitation light, you can see that they primarily
% escape in downward-pointing directions (in transmission), while the far
% field distribution of fluorescence light indicates that fluorescence
% mostly escapes on the ends of the cylinder, with a small amount of light
% coming out of the cylinder in the minus z direction (upwards).
%
% Note that MCmatlab distinguishes between "escaped" photons and "killed"
% photons. An "escaping" photon is one that hits the top cuboid boundary
% (if boundaryType == 2) or any cuboid boundary (if boundaryType == 1)
% where the medium has refractive index 1. A "killed" photon is one that
% strays too far from the main cuboid (5 times further than the cuboid
% dimensions).
%
% If boundaryType == 1, there are no "killed" photons since no photons can
% travel outside the cuboid, and the fraction of light absorbed in the
% cuboid plus the fraction of light escaping equals 1.
%% 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 = 101; % 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.simulationTimeRequested = .1; % [min] Time duration of the simulation
model.MC.farFieldRes = 50; % (Default: 0) If nonzero, photons that "escape" will have their energies tracked in a 2D angle distribution (theta,phi) array with theta and phi resolutions equal to this number. An "escaping" photon is one that hits the top cuboid boundary (if boundaryType == 2) or any cuboid boundary (if boundaryType == 1) where the medium has refractive index 1.
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 = 450; % [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 = .015; % [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 = 15/180*pi; % [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; % [cm] y position of focus
model.MC.lightSource.zFocus = 0.03; % [cm] z position of focus
model.MC.lightSource.theta = pi/6; % [rad] Polar angle of beam center axis
model.MC.lightSource.phi = -pi/2; % [rad] Azimuthal angle of beam center axis
model = runMonteCarlo(model);
model = plot(model,'MC');
figure(9); % Focus on the far field
%% Fluorescence Monte Carlo
model.FMC.simulationTimeRequested = 0.1; % [min] Time duration of the simulation
model.FMC.farFieldRes = 50; % (Default: 0) If nonzero, photons that "escape" will have their energies tracked in a 2D angle distribution (theta,phi) array with theta and phi resolutions equal to this number. An "escaping" photon is one that hits the top cuboid boundary (if boundaryType == 2) or any cuboid boundary (if boundaryType == 1) where the medium has refractive index 1.
model.FMC.matchedInterfaces = true; % Assumes all refractive indices are the same
model.FMC.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.FMC.wavelength = 550; % [nm] Fluorescence wavelength, used for determination of optical properties for fluorescence light
model = runMonteCarlo(model,'fluorescence');
model = plot(model,'FMC');
figure(19); % Focus on the far field
%% Geometry function(s) (see readme for details)
function M = geometryDefinition(X,Y,Z,parameters)
cylinderradius = 0.0100;
M = 1*ones(size(X)); % fill background with fluorescence absorber
M(Y.^2 + (Z - 3*cylinderradius).^2 < cylinderradius^2) = 2; % fluorescer
end
%% Media Properties function (see readme for details)
function mediaProperties = mediaPropertiesFunc(parameters)
mediaProperties = MCmatlab.mediumProperties;
j=1;
mediaProperties(j).name = 'fluorescence absorber';
mediaProperties(j).mua = @func1; % [cm^-1]
function mua = func1(wavelength)
if(wavelength<500)
mua = 1; % [cm^-1]
else
mua = 100; % [cm^-1]
end
end
mediaProperties(j).mus = 100; % [cm^-1]
mediaProperties(j).g = 0.9;
j=2;
mediaProperties(j).name = 'fluorescer';
mediaProperties(j).mua = @func2; % [cm^-1]
function mua = func2(wavelength)
if(wavelength<500)
mua = 100; % [cm^-1]
else
mua = 1; % [cm^-1]
end
end
mediaProperties(j).mus = 100; % [cm^-1]
mediaProperties(j).g = 0.9;
mediaProperties(j).QY = 0.4; % Fluorescence quantum yield
end