/mc_mainloop_rev4.m

http://github.com/gallamine/Photonator
Objective C | 356 lines | 313 code | 43 blank | 0 comment | 26 complexity | f777f29354897914b69b3acf78e3ab61 MD5 | raw file
```  1% Rev 4 - switch to split mover/receiver, fix direction vector problem,
2% lookup table/interpolation
3
4% set up an array for the photons,
5% x, y, z, mu_x, mu_y, mu_z, weight, received
6% 0, 0, 0, 0   , 0   , 1   , 1     , 1  - initial values
7% 1, 2, 3, 4   , 5   , 6   , 7     , 8  - Position
8% photon(:,1) == X POSITION
9% photon(:,2) == Y POSITION
10% photon(:,3) == Z POSITION
11% photon(:,4) == ux
12% photon(:,5) == uy
13% photon(:,6) == uz
14% photon(:,7) == WEIGHT
15% photon(:,8) == RECEIVED? 1 = No, 0 = Yes
16
17% index from origin, 0 theta (angle between x and z), 0 phi (angle between x and y) along x-axis
18
19clear all
20clc
21
22% Program run constants
23num_photons = 5e5;
24albedo = 0.95;          % Albedo of Maalox (ranges from 0.8 to 0.95)
25c = 2.35;
29
30rec_pos = [0,0; 0,0; 0,0];   % receivers at the same position (on the x/y plane)
31rec_aperture = [0.0079, 0.0508, 0.0508]; % rec_aperture = ones(num_rx,1).*0.0508;  % 0.0508 m = 2 inches
32rec_fov = [2.27, 0.0785, 0.1745];     % 0.314159 = 18 deg FOV, 2.27 = 130 deg, 0.0785 = 5 deg
33
34
35
36photon = zeros(num_photons,8);
37photon(:,7) = ones(num_photons,1);  % set weights to one
38photon(:,8) = ones(num_photons,1);  % set all active
39
40totaldist = zeros(num_photons,1);   % record total distance traveled by each photon
41rec_dist = zeros(num_photons,1);    % total distance photon traveld at point of reception
42rec_loc = zeros(num_photons,2);     % location of the received photon on the z,y rec. plane
43total_rec_packets = 0;              % Total number of packets to cross detector
44total_rec_power = 0;                % Total power to cross detector
45
46
47% c = 0.5;                            % attenuation coefficient in m^-1
48% a = 0.07;                           % absorption coefficient in m^-1
49% % harbor water
50% c = 2.190;
51% a = 0.366;
52% [cdf_scatter,angle] = generate_scatter('calc','harbor');
53
54% % coastal water
55% c = 0.22 + 0.179;
56% a = 0.179;
57% [cdf_scatter,angle] = generate_scatter('calc','coastal');
58% Clear water
59% c = 0.0374 + 0.114;
60% a = 0.114;
61% [cdf_scatter,angle] = generate_scatter('calc','clear');
62
63% Calculated values
64b = c * albedo;
65a = c - b;
66scattering_events = ceil((c-a)*receiver_z*5)     %five times the scattering attenuation length
67if scattering_events < 10
68      scattering_events = 10;
69end
70
71num_rx = length(rec_pos);
72
73inv_c = 1/c;
74inv_b = 1/(c-a);
75
76% Constants and magic numbers
77max_uz = 0.99999;
78n_water = 1.33;                     % index of refraction of water
79
80[cdf_scatter,angle] = generate_scatter('measured','maalox_alan');
81
82i = 0;
83
84h = waitbar(0.0,'Please wait...','CreateCancelBtn','stop=true; delete(h); clear h');
85set(h,'Name','optional window name');
86
87scattering_length = 1;
88
89% X position of the receiver (in meters)
90% receiver_z = 3*inv_c;                     % X position of the receiver (in meters)
91
93% a = c.*(1-albedo);
94
95% Photons at origin
96photon(:,1) = zeros(num_photons,1);        % 0
97photon(:,2) = zeros(num_photons,1);        % 0
98photon(:,3) = zeros(num_photons,1);        % 0
99% point down z-axis
100photon(:,4) = zeros(num_photons,1);        % 0
101photon(:,5) = zeros(num_photons,1);        % 0
102photon(:,6) = ones(num_photons,1);         % 1
103
104r_avg = 0;
105
106tic;
107
108for j = 1:scattering_events
109    rand_array = rand(num_photons,3);    % generate a matrix for each photon with rand propogation, roll, and pitch
110    rand_array(:,3) = rand_array(:,3).*2.*pi;   %% Uniformly distributed over 0 -2Pi
111    waitbar(j/scattering_events,h,['Scattering event ' num2str(j)]);
112
113    % iterate over every single photon to calculate new position and
114    % whether it was received or not.
115    for i = 1:num_photons
116
117        % if photon hasn't been received
118        if (photon(i,8) == 1)
119
120            r = -inv_c*log(rand_array(i,1));     % randomly generate optical path length
121
122
123             %Generate scattering angle from beam spread function
124            minIndex = 2;
125            maxIndex = length(cdf_scatter);
126            midIndex = minIndex + ceil((maxIndex - minIndex)/2);
127
128            while rand_array(i,2) ~= cdf_scatter(midIndex) && maxIndex > minIndex
129
130                midIndex = minIndex + ceil((maxIndex - minIndex)/2);
131                if rand_array(i,2) > cdf_scatter(midIndex)
132                    minIndex = midIndex + 1;
133                elseif rand_array(i,2) < cdf_scatter(midIndex)
134                    maxIndex = midIndex - 1;
135                end
136
137            end
138            midIndex = minIndex + ceil((maxIndex - minIndex)/2);
139            k = midIndex;
140
141            theta = angle(k-1) + (rand_array(i,2) - cdf_scatter(k-1))*(angle(k) - angle(k-1))/(cdf_scatter(k) - cdf_scatter(k-1));   % Linear interpolation from "angle" vector
142            phi = rand_array(i,3);  % Phi is uniformly distributed over 0-2pi
143
144
145            % find new position increment based on PREVIOUS direction vectors (ux,uy,uz). This
146            % takes care of the initial condition problem where photons
147            % were scattered immediately on the first iteration.
148            x_step = r*photon(i,4);  % x step size
149            y_step = r*photon(i,5);  % y step size
150            z_step = r*photon(i,6);  % z step size
151
152
153            % if the photon moves past the receiver plane
154            if ((photon(i,3) + z_step) >= receiver_z)
155
156%                 current_theta = acos(photon(i,6));      % cos^-1(uz)
157%                 if photon(i,6) ~= 1
158%                     if photon(i,5) >= 0
159%                         current_phi = acos(photon(i,4)/sqrt(1-photon(i,6)^2));  % cos^-1(ux/sqrt(1-uz^2))
160%                     else
161%                         current_phi = 2*pi - acos(photon(i,4)/sqrt(1-photon(i,6)^2));  % cos^-1(ux/sqrt(1-uz^2))
162%                     end
163%                 else
164%                     current_phi = rand()*2.*pi;
165%                 end
166%                 if ~isreal(current_phi)
167%                     format long
168%                     current_phi
169%                     photon(i,:)
170%                     disp('Invalid phi');
171%                     current_phi = real(current_phi);
172%                 end
173
174                % FIX THESE EQUATIONS!!
175
176                % z distance to receiver plane
177                z_dist_rec_intersection = receiver_z - photon(i,3);
178
179                if photon(i,6) ~= 0
180                    % y distance to receiver plane
181                    y_dist_rec_intersection = z_dist_rec_intersection*photon(i,5)/photon(i,6);      % z * uy/uz = z*tan(theta)*sin(phi)
182                    % x distance to receiver plane
183                    x_dist_rec_intersection = z_dist_rec_intersection*photon(i,4)/photon(i,6);      % z*ux/uz = z*tan(theta)*cos(phi)
184                else
185                    disp('how can the photon cross the receiver plane when it"s pointing straight up???');
186                end
187
188
189                % euclidian distance to the reciever plane
190                %dist_to_rec = sqrt((x_dist_rec_intersection)^2 + (y_dist_rec_intersection)^2 + (z_dist_rec_intersection)^2);
191                dist_to_rec = z_dist_rec_intersection / photon(i,6);    % z / mu_z
192
193
194                rec_loc(i,1) = photon(i,1) + x_dist_rec_intersection;   % x-axis location of reception
195                rec_loc(i,2) = photon(i,2) + y_dist_rec_intersection;    % y-axis location of reception
196                rec_loc(i,3) = acos(photon(i,6));                             % incident angle
197                rec_loc(i,4) = photon(i,4);                             % for statistics, should be uniform
198
199                total_rec_packets = total_rec_packets + 1;
200%                     total_rec_power = total_rec_power + photon(i,7)*exp(-dist_to_rec*a);
201                total_rec_power = total_rec_power + photon(i,7);
202                rec_dist(i) = totaldist(i)+ dist_to_rec;
203                photon(i,8) = 0;
204
205                totaldist(i) = totaldist(i) + dist_to_rec;
206
207
208
209            else        % if the photon didn't move into the detector, reduce its power & move it
210                % move to new x position
211                photon(i,1) = photon(i,1) + x_step;
212                % move to new y position
213                photon(i,2) = photon(i,2) + y_step;
214                % move to new z position
215                photon(i,3) = photon(i,3) + z_step;
216
217                if abs(photon(i,6)) > max_uz         % if uz ~ 1
218                    %photon(i,:)
219                    photon(i,4) = sin(theta)*cos(phi);
220                    photon(i,5) = sin(theta)*sin(phi);
221                    photon(i,6) = sign(photon(i,6))*cos(theta);
222                    %disp('mu_z near 1')
223                else
224                    sqrt_uz = sqrt(1 - photon(i,6)^2);
225                    old_ux = photon(i,4);
226                    old_uy = photon(i,5);
227                    old_uz = photon(i,6);
228                    photon(i,4) = (sin(theta)/sqrt_uz)*(old_ux*old_uz*cos(phi) - old_uy*sin(phi)) + old_ux*cos(theta);   % ux
229                    photon(i,5) = (sin(theta)/sqrt_uz)*(old_uy*old_uz*cos(phi) + old_ux*sin(phi)) + old_uy*cos(theta);   % uy
230                    photon(i,6) = (-sin(theta)*cos(phi))*sqrt_uz + old_uz*cos(theta);                                    % uz
231                    %disp('mu_z less than 1')
232                end
233
234                % update the total distance the photon has traveled
235                totaldist(i) = totaldist(i) + r;
236            end
237        end
238    end
239    %r_avg = r_avg/num_photons;
240%     figure(1)
241%     subplot(3,3,j)
242%     scatter(photon(:,1),photon(:,2),50,log10(photon(:,7)),'.')
243%     xlim([-15 50]);
244%     ylim([-25 25]);
245% %     xlabel('x-axis (m)')
246% %     xlabel('y-axis (m)')
247end
248
249
250
251
252
253rec_loc_final = ones(total_rec_packets,4);
254j = 1;
255for i = 1:num_photons
256    if (photon(i,8) == 0)
257       rec_loc_final(j,:) = rec_loc(i,:);
258       j = j +1;
259    end
260end
261
262j = 1;
263total_rec_dist = zeros(total_rec_packets,1);
264total_rec_power = 0;
265for i = 1:num_photons
266
267   if (photon(i,8) == 0)
268       total_rec_dist(j) = rec_dist(i);
269       total_rec_power = total_rec_power + exp(-a*rec_dist(i));
270       j = j + 1;
271   end
272end
273
274
275[power,ph_cnt,angle_mean,angle_var,power_mean,power_var] = mc_rec_r3(a,rec_loc_final,total_rec_dist,rec_pos,rec_aperture,rec_fov);
276
277power/num_photons               % Total received power / total transmitted photon packets
278sd_power = sqrt(power_var).*ph_cnt/num_photons;     % sd_power = sd_photon*rx_photons/total_photons.
279
280
281total_time = toc;
282minutes_to_run = total_time/60
283
284% figure(8)                   % plot 3D histogram of photons on RX plane
285% hist3([rec_loc_final(:,1),rec_loc_final(:,2)],[100 100])
286% set(gcf,'renderer','opengl');
287% set(get(gca,'child'),'FaceColor','interp','CDataMode','auto');
288
289% j = 1;
290% scattered_times = zeros(total_rec_packets,1);
291% for i = 1:num_photons
292%
293%    if (photon(i,7) == 0)
294%        scattered_times(j) = photon(i,8);
295%        j = j + 1;
296%    end
297% end
298
299total_rec_power
300total_rec_packets
301% hist(scattered_times,[0:10]);
302
303close(h)
304
305% figure(10)                              % Plot histogram of time-of-arrival vs. power
306% [N,X] = hist(total_rec_dist,100);
307% pwr_rx_hist = N.*exp(-a.*X);
308% stem(X/(3e8/n_water),pwr_rx_hist/total_rec_power)
309% xlabel('Time of arrival (sec)')
310% ylabel('Percentage of power')
311
312% distance_delta = max(X) - min(X);
313% time_delta = distance_delta/(3e8/n_water);
314% T = mean(X(2:end) - X(1:end-1))/(3e8/n_water); % Effective "sampling rate" of the histogram
315% bandwidth = 1/T;                                % Normalized frequency in Hz
316
317%  figure(7)
318% freqz(N/total_rec_packets,[1],512,bandwidth)      %% plot a frequency response from the histogram data
319
320% figure(9)                           %% Scatter plot of photons on RX plane
321% scatter(rec_loc(:,1),rec_loc(:,2))
322
323% figure(1)
324% scatter(photon(:,1),photon(:,2),50,log10(photon(:,7)),'.')
325% xlim([-20 150]);
326% ylim([-120 120]);
327% xlabel('x-axis (m)')
328% xlabel('y-axis (m)')
329%
330% figure(4)
331% scatter3(photon(:,1),photon(:,2),photon(:,3),50,log10(photon(:,7)),'.')
332
333% % Draw the box representing the receiver
338% xlabel('x-axis (m)')
339% ylabel('y-axis (m)')
340% zlabel('z-axis (m)')
341
342% figure(3);
343% hist(photon(:,1),max(photon(:,1)));
344
345
346% figure(2)
347% hist(totaldist,20);
348% figure(3)
349% hist(photon(:,4),20)
350findfigs
351
352sprintf('Simulation on DATE with %d photons, %d scattering events.', num_photons, scattering_events)
353sprintf('C = %d (1/m), A = %d (1/m).',c,a)