/mc_func_Berrocal.m
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