/mc_rec_r4.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 <- 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