PageRenderTime 69ms CodeModel.GetById 14ms app.highlight 52ms RepoModel.GetById 1ms app.codeStats 0ms

/mc_func_Berrocal.m

http://github.com/gallamine/Photonator
Objective C | 249 lines | 219 code | 30 blank | 0 comment | 26 complexity | ea7e5e033ce090ef6a793b1193e4dfeb MD5 | raw file
  1% set up an array for the photons, 
  2% x, y, z, mu_x, mu_y, mu_z, weight, received
  3% 0, 0, 0, 0   , 0   , 1   , 1     , 1  - initial values
  4% 1, 2, 3, 4   , 5   , 6   , 7     , 8  - Position 
  5% photon(:,1) == X POSITION
  6% photon(:,2) == Y POSITION
  7% photon(:,3) == Z POSITION
  8% photon(:,4) == ux
  9% photon(:,5) == uy
 10% photon(:,6) == uz
 11% photon(:,7) == WEIGHT
 12% photon(:,8) == RECEIVED? 1 = No, 0 = Yes (detected), -1 = terminated
 13
 14% index from origin, 0 theta (angle between x and z), 0 phi (angle between
 15% x and y) along x-axis
 16
 17% V5 - change scattering equations. Change direction of propogation to
 18% z-axis.
 19% V6 - Lots of changes:
 20%   - Added size limits to the simulation. Photons terminate when
 21%   encountering boundaris
 22%   - Added unlimited scattering. Scattering is dependant on minimum weight
 23%   - added rouletting to remove photons of low weight
 24%   - added window effects to the receiver
 25%   - changed events to be based on 'c', not 'b'
 26%   - changed how weights are calculated/updated
 27
 28function [total_time,total_rec_power,total_rec_packets,rec_loc_final,total_rec_dist,rec_weights] = ...
 29    mc_func_Berrocal(num_photons,scattering_events,c,a,receiver_z,...
 30    cdf_scatter,angle,init_angle,init_angle2,diverg)
 31
 32useLimits = 'true';
 33
 34if strcmp(useLimits,'true')
 35    xLimMax = 0.005000;
 36    xLimMin = -0.005;
 37    yLimMax = 0.005;
 38    yLimMin = -0.005;
 39    zLimMax = 0.01;
 40    zLimMin = -0.005;
 41end
 42
 43photon = zeros(num_photons,8);
 44photon(:,7) = ones(num_photons,1);  % set weights to one
 45photon(:,8) = ones(num_photons,1);  % set all active
 46totaldist = zeros(num_photons,1);   % record total distance traveled by each photon
 47rec_dist = zeros(num_photons,1);    % total distance photon traveld at point of reception
 48rec_loc = zeros(num_photons,4);     % location of the received photon on the z,y rec. plane
 49total_rec_packets = 0;              % Total number of packets to cross detector
 50total_rec_power = 0;                % Total power to cross detector
 51
 52prob_of_survival = ((c-a)/c);       % b/c
 53rouletteConst = 10;                 % 1/rouletteConst determines the probability of a low power photon surviving
 54rouletteConst_inv = 1/rouletteConst;
 55inv_c = 1/c;
 56inv_b = 1/(c-a);
 57
 58min_power = 1e-4;                   % minimum power value for photons before they are terminated by rouletting
 59
 60max_uz = 0.99999;
 61
 62tic;
 63% 
 64% theta = diverg.*rand(1,num_photons);    % uniform rand. dist over divergence ang.
 65% phi = (2*pi).*rand(1,num_photons);      % uniform rand. dist over azmuthal angles
 66% 
 67% beta = init_angle;   % direction the transmitter is pointing (zenith)
 68% alpha = init_angle2;     % direction the transmitter is pointing (azmuth)
 69
 70
 71[x,y] = startingDistribution(num_photons);         % Assuming a square end size
 72
 73%point down z-axis
 74photon(:,1) = x.*xLimMax;                              % x position
 75photon(:,2) = y.*yLimMax;                           % y position
 76photon(:,6) = ones(num_photons,1);         % z - 1
 77
 78% photon(:,1) = xLimMax;                              % x position
 79% photon(:,2) = y.*yLimMax;                           % y position
 80% photon(:,3) = x.*xLimMax;                           % z position
 81% photon(:,4) = -1.*ones(num_photons,1);         % x - 1
 82
 83photonsRemaining = num_photons;             % count down for each photon received/terminated
 84
 85clear theta phi  z beta alpha
 86
 87% for j = 1:scattering_events
 88while photonsRemaining > 0                      % which is faster? create random values on the fly??  <- check
 89    rand_array = rand(num_photons,3);           % generate a matrix for each photon with rand propogation, roll, and pitch
 90    rand_array(:,3) = rand_array(:,3).*2.*pi;   % Uniformly distributed over (0,2Pi)                                                                             
 91    
 92    % iterate over every single photon to calculate new position and
 93    % whether it was received or not.
 94    for i = 1:num_photons
 95        
 96        % if photon hasn't been received
 97        if (photon(i,8) == 1)
 98        
 99             r = -inv_c*log(rand_array(i,1));     % randomly generate optical path length      
100
101            %Generate scattering angle from beam spread function
102            minIndex = 2;
103            maxIndex = length(cdf_scatter);
104            midIndex = minIndex + ceil((maxIndex - minIndex)/2);
105
106            while rand_array(i,2) ~= cdf_scatter(midIndex) && maxIndex >= minIndex        % changed > to >=      
107    
108                midIndex = minIndex + ceil((maxIndex - minIndex)/2);
109                if rand_array(i,2) > cdf_scatter(midIndex)
110                    minIndex = midIndex + 1;
111                elseif rand_array(i,2) < cdf_scatter(midIndex)
112                    maxIndex = midIndex - 1;
113                end
114
115            end
116            midIndex = minIndex + ceil((maxIndex - minIndex)/2);
117            k = midIndex;
118
119            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
120            phi = rand_array(i,3);  % Phi is uniformly distributed over 0-2pi
121
122
123            % find new position increment based on PREVIOUS direction vectors (ux,uy,uz). This
124            % takes care of the initial condition problem where photons
125            % were scattered immediately on the first iteration.
126            x_step = r*photon(i,4);  % x step size
127            y_step = r*photon(i,5);  % y step size
128            z_step = r*photon(i,6);  % z step size
129                
130
131            % if the photon moves past the receiver plane
132            if ((photon(i,3) + z_step) >= receiver_z)
133          
134                if photon(i,6) ~= 0                         % If the photon has a z-component, mu_z != 0
135                    % z distance to receiver plane
136                    z_dist_rec_intersection = receiver_z - photon(i,3);     % Zrec - Zphoton
137                    % y distance to receiver plane
138                    y_dist_rec_intersection = z_dist_rec_intersection*photon(i,5)/photon(i,6);
139                    % x distance to receiver plane
140                    x_dist_rec_intersection = z_dist_rec_intersection*photon(i,4)/photon(i,6);
141                else
142                    disp('how can the photon cross the receiver plane when it"s pointing straight up??? Z distance step:');
143                    disp(z_step);
144                    disp('Z position:');
145                    disp(photon(i,3));
146                end
147
148                % euclidian distance to the reciever plane
149                dist_to_rec = z_dist_rec_intersection / photon(i,6);    % z / mu_z
150
151                rec_loc(i,1) = photon(i,1) + x_dist_rec_intersection;   % x-axis location of reception
152                rec_loc(i,2) = photon(i,2) + y_dist_rec_intersection;    % y-axis location of reception
153                rec_loc(i,3) = photon(i,4);                             % incident angle (mu_x)
154                rec_loc(i,4) = photon(i,5);                             % for statistics, should be uniform (mu_y)
155                rec_loc(i,5) = photon(i,6);                             % for statistics, should be uniform (mu_z)
156
157                total_rec_packets = total_rec_packets + 1;
158                total_rec_power = total_rec_power + photon(i,7);        % total power at the receiver plane (sum of received photons)
159                
160                rec_dist(i) = totaldist(i)+ dist_to_rec;                % individual photon's distance traveled.
161                photon(i,8) = 0;                                        % mark photon as received
162                photonsRemaining = photonsRemaining - 1;                % decrement number of outstanding photons
163                
164                % update the total distance the photon has traveled 
165                totaldist(i) = totaldist(i) + dist_to_rec;
166
167            else % if the photon didn't move into the detector, reduce its power & move it
168                
169                photon(i,1) = photon(i,1) + x_step;         % move to new x position                            
170                photon(i,2) = photon(i,2) + y_step;         % move to new y position                        
171                photon(i,3) = photon(i,3) + z_step;         % move to new z position
172                
173                if ((photon(i,1) > xLimMax) || (photon(i,1) < xLimMin) || (photon(i,2) > yLimMax) || (photon(i,2) < yLimMin) || (photon(i,3) < zLimMin))
174                    photon(i,8) = -1;                           % mark as terminated
175                    photonsRemaining = photonsRemaining - 1;    % decrement outstanding photons
176                else 	% if the photon is still in the boundaries                      
177               
178                    photon(i,7) = photon(i,7)*prob_of_survival;     % reduce weight
179                    if  photon(i,7) < min_power
180                        if rand() > (rouletteConst_inv)                 % if the photon is randomly chosen to be terminated ...
181                            photon(i,8) = -1;                               % -1 indicates a photon was terminated, but not received
182                            photonsRemaining = photonsRemaining - 1;        % decrement outstanding photons
183                        else                                            % ... otherwise the photon gets the energey of terminated photons
184                            photon(i,7) = photon(i,7)*rouletteConst;    % shift power of terminated photons to new photon
185                        end
186                    end
187                    
188                    if abs(photon(i,6)) > max_uz         % if uz ~ 1
189                        photon(i,4) = sin(theta)*cos(phi);
190                        photon(i,5) = sin(theta)*sin(phi);
191                        photon(i,6) = sign(photon(i,6))*cos(theta);
192                        %disp('mu_z near 1')
193                    else
194                        sqrt_uz = sqrt(1 - photon(i,6)^2);
195                        old_ux = photon(i,4);
196                        old_uy = photon(i,5);
197                        old_uz = photon(i,6);
198                        photon(i,4) = (sin(theta)/sqrt_uz)*(old_ux*old_uz*cos(phi) - old_uy*sin(phi)) + old_ux*cos(theta);   % ux
199                        photon(i,5) = (sin(theta)/sqrt_uz)*(old_uy*old_uz*cos(phi) + old_ux*sin(phi)) + old_uy*cos(theta);   % uy
200                        photon(i,6) = (-sin(theta)*cos(phi))*sqrt_uz + old_uz*cos(theta);                                    % uz
201                        %disp('mu_z less than 1')
202                    end
203                end
204                % update the total distance the photon has traveled 
205                totaldist(i) = totaldist(i) + r;
206
207            end
208        end
209    end
210end
211
212total_time = toc;
213
214% Do this after the loop above, so we can allocate the array once (instead
215% of growing dynamically as we receive individual photons - SLOW)
216rec_loc_final = ones(total_rec_packets,5);
217j = 1;
218for i = 1:num_photons                           % iterate over all photons
219    if (photon(i,8) == 0)                       % if the photon was received
220       rec_loc_final(j,:) = rec_loc(i,:);       % record the receive location and angles
221       j = j +1;                                % increment the number of received photons
222    end
223end
224
225if ((j-1) ~= total_rec_packets)
226    disp('Error! Total received photons doesnt equal j iterator. ');
227    disp(sprintf('j = %d and total_rec_packets = %d',j, total_rec_packets));
228end
229
230j = 1;
231total_rec_dist = zeros(total_rec_packets,1);
232rec_weights = zeros(total_rec_packets,1);
233total_rec_power = 0;
234for i = 1:num_photons    
235   if (photon(i,8) == 0)
236       total_rec_dist(j) = rec_dist(i);         % store the distance traveled by each photon
237       rec_weights(j) = photon(i,7);            % store the weights of received photons
238       j = j + 1;       
239   end
240end
241
242if ((j-1) ~= total_rec_packets)
243    disp('Error! Total received distances doesnt equal j iterator');
244    disp(sprintf('j = %d and total_rec_packets = %d',j, total_rec_packets));
245end
246
247
248
249