/mc_rec_r4.m
Objective C | 172 lines | 138 code | 34 blank | 0 comment | 30 complexity | d76a5981c58da01c749e3d5c5ad06f3f MD5 | raw file
1% VERSION LOG 2% Rev 1 - first creation, record power and count for each receiver 3% Rev 2 - add statistics, mean and variance for power, angle(s) 4% Rev 3 - changed receiver plane to x/y instead of y/z - photons move along 5% the z-axis now 6% Rev 3b/4 - 5/3/11 - added window affects, changed power calculations <- SIMPLE BASIC, VERIFIED VERSION 7 8function [power,ph_cnt,angle_mean,angle_var,dist_mean,dist_var,weight_mean,weight_var,reflected,distances,angles,weights] = mc_rec_r4(a,rec_loc_final,total_rec_dist,rec_weights,rec_pos,rec_aperture,rec_fov,numTxPhotons) 9 10sizeRecPos = size(rec_pos); 11num_rx = sizeRecPos(1); % total number of receivers 12 13num_photons = length(rec_loc_final); % total number of received photons 14 15ph_cnt = zeros(1,num_rx); 16power = zeros(1,num_rx); 17 18angle_mean = zeros(1,num_rx); % mean of the incident angle for each Rx 19angle_var = zeros(1,num_rx); % variance of incident angle for each Rx 20dist_mean = zeros(1,num_rx); % Mean distance traveled for each Rx 21 22dist_var = zeros(1,num_rx); % Variance of dist traveled for each Rx 23weight_mean = zeros(1,num_rx); % Mean weight for each Rx 24weight_var = zeros(1,num_rx); % Variance of weight for each Rx 25 26distances = zeros(1,num_photons); % Distances from photon to 0,0 point - Initailize to full size, crop at end 27angles = zeros(1,num_photons); % Angle of received photon 28weights = zeros(1,num_photons); % Weight of received photon 29 30dWindow = 0.00635; % 0.25 inches 31dAir = 0.127; % 5 inches 32nWater = 1.33; 33nWindow = 1.585; % polycarbonate 34nAir = 1; 35reflected = 0; % number of photons past critical angle of window 36nWaterWindow = nWater/nWindow; 37nWindowAir = nWindow/nAir; 38 39critAngSine = nAir/nWater; % sine T1 = nair/nwater sine T3 40 41if (num_rx > 1) 42 disp('You need to change how the distance array is stored!'); 43end 44 45for i = 1:num_photons % iterate over every photon on receiver plane 46 ph_x = rec_loc_final(i,1); 47 ph_y = rec_loc_final(i,2); 48 mu_x = rec_loc_final(i,3); 49 mu_y = rec_loc_final(i,4); 50 mu_z = rec_loc_final(i,5); 51 52 % Find new position from window 53 sinT = sqrt(1 - mu_z^2); 54 55% if sinT >= critAngSine % if the photon's incident angle > critical angle, then skip this photon 56% reflected = reflected + 1; 57% continue; 58% end 59% 60% sinT1 = nWaterWindow*sinT; % 1.33/1.585 * sqrt(1 - uz^2) 61% mu_z1 = sqrt(1-sinT1^2); % this can be simplified ... 62% mu_x1 = nWaterWindow * mu_x; 63% mu_y1 = nWaterWindow * mu_y; 64% delY1 = dWindow * mu_y1/mu_z1; 65% delX1 = dWindow * mu_x1/mu_z1; 66% % find new position from air 67% sinT2 = nWater*sinT; % this can be simplified ... 68% mu_z2 = sqrt(1 - sinT2^2); 69% mu_x2 = nWater * mu_x; 70% mu_y2 = nWater * mu_y; 71% delY2 = dAir * mu_y2/mu_z2; 72% delX2 = dAir * mu_x2/mu_z2; 73% 74% ph_x = ph_x + delX1 + delX2; 75% ph_y = ph_y + delY1 + delY2; 76 77 for j = 1:num_rx % iterate over every receiver 78 rx_x = rec_pos(j,1); 79 rx_y = rec_pos(j,2); 80 radius = rec_aperture(j)/2; % 1/2 diameter of receiver 81 angle_sqravg = 0; 82 dist_sqravg = 0; 83 84 cos_rec_fov = cos(rec_fov(j)/2); % cos(fov/2) to compare with photon's incident angle 85 86 distance = sqrt((ph_x-rx_x)^2 + (ph_y-rx_y)^2); % Euclidian distance to receiver center 87 88 if ((distance <= radius) && (mu_z >= cos_rec_fov)) % Photon received 89 power(j) = power(j) + rec_weights(i); % Adjust power of received photons 90 ph_cnt(j) = ph_cnt(j) + 1; % Increment number of received photons 91 ph_dist = total_rec_dist(i); 92 93 angle_delta = rec_loc_final(i,5) - angle_mean(j); 94 dist_delta = ph_dist - dist_mean(j); 95 weight_delta = rec_weights(i) - weight_mean(j); 96 97% if j == 1 98% weights(ph_cnt(j)) = rec_weights(i); 99% angles(ph_cnt(j)) = rec_loc_final(i,5); 100% distances(ph_cnt(j)) = ph_dist; 101% end 102 103 % Record stats on incident angle 104 angle_mean(j) = angle_mean(j) + (angle_delta) / ph_cnt(j); % E'[uz] = E[uz] + (uz - E[uz]) / N+1 - See Knuth's Art of Comp. Programming & http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance 105 % Record stats on distance traveled 106 dist_mean(j) = dist_mean(j) + (dist_delta) / ph_cnt(j); 107 % Record stats on photon weight 108 weight_mean(j) = weight_mean(j) + (weight_delta)/ph_cnt(j); % E'[D] = E[D] + (D - E[D]) / N+1 109 110 111 angle_var(j) = angle_var(j) + angle_delta*(rec_loc_final(i,5) - angle_mean(j)); % Var[uz] = 1/N-1 *(Var[uz] + (uz - E[uz]_N-1)(uz - E[uz]_N)) 112 dist_var(j) = dist_var(j) + dist_delta*(ph_dist - dist_mean(j)); 113 weight_var(j) = weight_var(j) + weight_delta*(rec_weights(i) - weight_mean(j)); 114 115 if (num_rx == 1) 116 distances(ph_cnt(j)) = distance; 117 angles(ph_cnt(j)) = mu_z; 118 weights(ph_cnt(j)) = rec_weights(i); 119 end 120 end 121 end 122end 123 124 125 126for j = 1:num_rx 127 128 if (num_rx == 1) 129 distances = distances(1:ph_cnt(j)); 130 angles = angles(1:ph_cnt(j)); 131 weights = weights(1:ph_cnt(j)); 132 end 133 134 weight_var(j) = (weight_var(j) + weight_mean(j).^2*(ph_cnt(j).*(numTxPhotons - ph_cnt(j))/numTxPhotons)) / (numTxPhotons - 1); 135 weight_mean(j) = ph_cnt(j)*weight_mean(j) / numTxPhotons; 136 137 if weight_mean(j) ~= power(j)/numTxPhotons 138 disp('Problem counting up the weight mean'); 139 weight_mean(j) - power(j)/numTxPhotons 140 weight_mean(j) 141 power(j) 142 end 143 144% if ph_cnt(j) > 1 145% angle_var(j) = (1./(ph_cnt(j)-1)).*angle_var(j); 146% dist_var(j) = (1./(ph_cnt(j)-1)).*dist_var(j); 147% end 148% 149% if isnan(angle_var(j)) 150% disp('angle_var(j) is NaN (in correction loop)'); 151% end 152% if isnan(dist_var(j)) 153% disp('dist_var(j) is NaN (in correction loop)'); 154% end 155% if isnan(weight_var(j)) 156% disp('weight_var(j) is NaN (in correction loop)'); 157% end 158% if isinf(angle_var(j)) 159% disp('angle_var(j) is Inf'); 160% end 161% if isinf(dist_var(j)) 162% disp('dist_var(j) is Inf'); 163% end 164% if isinf(weight_var(j)) 165% disp('weight_var(j) is Inf'); 166% end 167 168end 169 170 171 172%dist_mean.*ph_cnt - power