/mc_rec_r5.m
Objective C | 214 lines | 170 code | 44 blank | 0 comment | 35 complexity | 6fe170be6a630c929684a7435c8e6477 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 7% Rev 5 - 7/22/11 - added FOV function instead of hard limits! Added window 8% effects 9 10function [power,ph_cnt,angle_mean,angle_var,dist_mean,dist_var,weight_mean,weight_var,reflected,distances,angles,weights] = mc_rec_r5(a,rec_loc_final,total_rec_dist,rec_weights,rec_pos,rec_aperture,rec_fov,numTxPhotons) 11 12% hard coded receiver function, one receiver right now 13% Two-term Gaussian 14% Terms 15 16A1 = 0.7985; 17B1 = 0.0187; 18C1 = 0.03437; 19A2 = 0.7121; 20B2 = -0.02337; 21C2 = 0.03117; 22 23sizeRecPos = size(rec_pos); 24if sizeRecPos > 1 25 error('Receiver V5 only supports one receiver. HARD CODED!!!!'); 26end 27num_rx = sizeRecPos(1); % total number of receivers 28 29num_photons = length(rec_loc_final); % total number of received photons 30 31ph_cnt = zeros(1,num_rx); 32power = zeros(1,num_rx); 33 34angle_mean = zeros(1,num_rx); % mean of the incident angle for each Rx 35angle_var = zeros(1,num_rx); % variance of incident angle for each Rx 36dist_mean = zeros(1,num_rx); % Mean distance traveled for each Rx 37 38dist_var = zeros(1,num_rx); % Variance of dist traveled for each Rx 39weight_mean = zeros(1,num_rx); % Mean weight for each Rx 40weight_var = zeros(1,num_rx); % Variance of weight for each Rx 41 42distances = zeros(1,num_photons); % Distances from photon to 0,0 point - Initailize to full size, crop at end 43angles = zeros(1,num_photons); % Angle of received photon 44weights = zeros(1,num_photons); % Weight of received photon 45 46dWindow = 0.00635; % 0.25 inches 47dAir = 0.127; % 5 inches 48nWater = 1.33; 49nWindow = 1.585; % polycarbonate 50nAir = 1; 51reflected = 0; % number of photons past critical angle of window 52nWaterWindow = nWater/nWindow; 53nWindowAir = nWindow/nAir; 54 55critAngCos = sqrt(1 - (nAir/nWater)^2); % sine T1 = nair/nwater sine T3 56 57 58albedo = 0.95; 59scatterLimit = -1; 60weightLimit = albedo^scatterLimit; 61 62if scatterLimit ~= -1 63 disp('The receiver will only receive photon scattered a fixed amount of times!'); 64end 65 66if (num_rx > 1) 67 disp('You need to change how the distance array is stored!'); 68end 69 70for i = 1:num_photons % iterate over every photon on receiver plane 71 ph_x = rec_loc_final(i,1); 72 ph_y = rec_loc_final(i,2); 73 mu_x = rec_loc_final(i,3); 74 mu_y = rec_loc_final(i,4); 75 mu_z = rec_loc_final(i,5); 76 77 % Find new position from window 78 79 if mu_z <= critAngCos % if the photon's incident angle > critical angle, then skip this photon 80 reflected = reflected + 1; 81 continue; 82 end 83 84 % Reject photons that are single/multiple scattered (used for testing 85 % purposes and to figure out where the power comes from) 86 if (scatterLimit ~= -1) && (rec_weights(i) < weightLimit) 87 continue; 88 end 89 90 % Calculate fresnel loss 91 cosExitAng = sqrt(1 - (nWater/nAir)^2*(1-mu_z^2)); 92 rp = (mu_z - nWater*cosExitAng)/(mu_z + nWater*cosExitAng); % Thanks to -> Small Monte Carlo by Scott Prahl (http://omlc.ogi.edu) 93 rs = (cosExitAng - nWater*mu_z)/(cosExitAng + nWater*mu_z); 94 R = (rp^2 + rs^2)/2; % unpolarized reflection coefficient 95 T = 1-R; % percent transmitted 96 97 rec_weights(i) = rec_weights(i)*T; 98 mu_z = cosExitAng; 99 100 % 101% sinT1 = nWaterWindow*sinT; % 1.33/1.585 * sqrt(1 - uz^2) 102% mu_z1 = sqrt(1-sinT1^2); % this can be simplified ... 103% mu_x1 = nWaterWindow * mu_x; 104% mu_y1 = nWaterWindow * mu_y; 105% delY1 = dWindow * mu_y1/mu_z1; 106% delX1 = dWindow * mu_x1/mu_z1; 107% % find new position from air 108% sinT2 = nWater*sinT; % this can be simplified ... 109% mu_z2 = sqrt(1 - sinT2^2); 110% mu_x2 = nWater * mu_x; 111% mu_y2 = nWater * mu_y; 112% delY2 = dAir * mu_y2/mu_z2; 113% delX2 = dAir * mu_x2/mu_z2; 114% 115% ph_x = ph_x + delX1 + delX2; 116% ph_y = ph_y + delY1 + delY2; 117 118 for j = 1:num_rx % iterate over every receiver 119 rx_x = rec_pos(j,1); 120 rx_y = rec_pos(j,2); 121 radius = rec_aperture(j)/2; % 1/2 diameter of receiver 122 123 cos_rec_fov = cos(rec_fov(j)/2); % cos(fov/2) to compare with photon's incident angle 124 125 distance = sqrt((ph_x-rx_x)^2 + (ph_y-rx_y)^2); % Euclidian distance to receiver center 126 127 if ((distance <= radius) && (mu_z >= cos_rec_fov)) % Photon received 128 arrivalAng = acos(mu_z); 129 fovWeight = A1*exp(-((arrivalAng-B1)/C1)^2) + A2*exp(-((arrivalAng-B2)/C2)^2); 130 power(j) = power(j) + rec_weights(i)*fovWeight; % Adjust power of received photons 131 132 ph_cnt(j) = ph_cnt(j) + 1; % Increment number of received photons 133 ph_dist = total_rec_dist(i); 134 135 angle_delta = rec_loc_final(i,5) - angle_mean(j); 136 dist_delta = ph_dist - dist_mean(j); 137 weight_delta = rec_weights(i)*fovWeight - weight_mean(j); 138 139% if j == 1 140% weights(ph_cnt(j)) = rec_weights(i)*fovWeight; 141% angles(ph_cnt(j)) = rec_loc_final(i,5); 142% distances(ph_cnt(j)) = ph_dist; 143% end 144 145 % Record stats on incident angle 146 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 147 % Record stats on distance traveled 148 dist_mean(j) = dist_mean(j) + (dist_delta) / ph_cnt(j); 149 % Record stats on photon weight 150 weight_mean(j) = weight_mean(j) + (weight_delta)/ph_cnt(j); % E'[D] = E[D] + (D - E[D]) / N+1 151 152 153 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)) 154 dist_var(j) = dist_var(j) + dist_delta*(ph_dist - dist_mean(j)); 155 weight_var(j) = weight_var(j) + weight_delta*(rec_weights(i)*fovWeight - weight_mean(j)); 156 157% if (num_rx == 1) 158% distances(ph_cnt(j)) = distance; 159% angles(ph_cnt(j)) = mu_z; 160% weights(ph_cnt(j)) = rec_weights(i)*fovWeight; 161% end 162 end 163 end 164end 165 166 167 168for j = 1:num_rx 169 170 if (num_rx == 1) 171 distances = distances(1:ph_cnt(j)); 172 angles = angles(1:ph_cnt(j)); 173 weights = weights(1:ph_cnt(j)); 174 end 175 176 weight_var(j) = (weight_var(j) + weight_mean(j).^2*(ph_cnt(j).*(numTxPhotons - ph_cnt(j))/numTxPhotons)) / (numTxPhotons - 1); 177 weight_mean(j) = ph_cnt(j)*weight_mean(j) / numTxPhotons; 178 179 if weight_mean(j) ~= power(j)/numTxPhotons 180 %disp('Problem counting up the weight mean'); 181 %weight_mean(j) - power(j)/numTxPhotons; 182 weight_mean(j); 183 power(j); 184 end 185 186% if ph_cnt(j) > 1 187% angle_var(j) = (1./(ph_cnt(j)-1)).*angle_var(j); 188% dist_var(j) = (1./(ph_cnt(j)-1)).*dist_var(j); 189% end 190% 191% if isnan(angle_var(j)) 192% disp('angle_var(j) is NaN (in correction loop)'); 193% end 194% if isnan(dist_var(j)) 195% disp('dist_var(j) is NaN (in correction loop)'); 196% end 197% if isnan(weight_var(j)) 198% disp('weight_var(j) is NaN (in correction loop)'); 199% end 200% if isinf(angle_var(j)) 201% disp('angle_var(j) is Inf'); 202% end 203% if isinf(dist_var(j)) 204% disp('dist_var(j) is Inf'); 205% end 206% if isinf(weight_var(j)) 207% disp('weight_var(j) is Inf'); 208% end 209 210end 211 212 213 214%dist_mean.*ph_cnt - power