/mc_rec_r5.m

  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