/mc_mainloop_rev4.m
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; 26receiver_z = 3.66; 27receiver_x = 0.20; % Y position of the receiver (in meters) 28receiver_y = 0; % Z position of the receiver(in meters) 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 92% c = scattering_length./receiver_z; 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 334% line([receiver_z receiver_z],[receiver_y_min receiver_y_min],[receiver_z_max receiver_z_min],'LineWidth',4) % | 335% line([receiver_z receiver_z],[receiver_y_min receiver_y_max],[receiver_z_max receiver_z_max],'LineWidth',4) % - 336% line([receiver_z receiver_z],[receiver_y_max receiver_y_max],[receiver_z_max receiver_z_min],'LineWidth',4) % | 337% line([receiver_z receiver_z],[receiver_y_min receiver_y_max],[receiver_z_min receiver_z_min],'LineWidth',4) % _ 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) 354sprintf('Receiver at %d, %d, %d (meters)',receiver_x, receiver_y, receiver_z) 355% sprintf('Travel distance delta %d (m). Time of arrival delta %d (sec)', distance_delta, time_delta) 356% sprintf('Time delta between histogram bins: %d (sec), %d (Hz)',T,bandwidth)