/mc_rec_r5.m

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