Rewrite of function with significant speed up

1 files changed, 29 insertions, 22 deletions

diff --git a/gbeam_propagation_froward_only.m b/gbeam_propagation_froward_only.m
index db8877f..93fe2c7 100644
--- a/gbeam_propagation_froward_only.m
+++ b/gbeam_propagation_froward_only.m
@@ -15,31 +15,38 @@ function q = gbeam_propagation_froward_only(x_pos, q_in, x_in,  optics_elements)
 	optics_elements=arrange_optics_along_x(optics_elements);
 
 	% Forward propagation to the right of x_in
-	Np=length(x_pos); % number of 'x' points
-	Nel=length(optics_elements) ;
 	q=0*x_pos; % q vector initialization
 	q_last_calc=q_in;
 	x_last_calc=x_in; % the furthest calculated point
-	for i=1:Np
-		x_pos_i=x_pos(i);
-		for k=1:length(optics_elements) 
-			% iterates through optics_elements to make sure 
-			% we take them in account for the beam propagation
-			el=optics_elements{k};
-			if ( (x_last_calc < el.x) && (el.x <= x_pos_i) )
-				abcd=abcd_free_space(el.x-x_last_calc);
-				q_last_calc=q_after_abcd(q_last_calc,abcd);
-				q_last_calc=q_after_abcd(q_last_calc,el.abcd);
-				x_last_calc=el.x;
-			end
-		end
-		if (x_pos_i > x_last_calc);
-			abcd=abcd_free_space(x_pos_i-x_last_calc);
-			q_last_calc=q_after_abcd(q_last_calc,abcd);
-			x_last_calc=x_pos_i;
-		end
-		q(i)=q_last_calc;
-	end
+    
+    
+    for k = 1:length(optics_elements)
+        el=optics_elements{k};
+        if el.x > x_last_calc
+            %Find all points to the left of the current optical element
+            index = (x_last_calc <= x_pos) & (x_pos < el.x);
+            if any(index)
+                x=x_pos(index);
+                q_tmp=q_after_free_space(q_last_calc,x-x_last_calc);
+                q(index)=q_tmp;
+                x_last_calc=x(end);
+                q_last_calc=q_tmp(end);
+            end
+            %propagate beam up to the lens
+            q_at_lens = q_after_free_space(q_last_calc,el.x - x_last_calc);
+            %Applying lens transformation
+            q_last_calc=q_after_abcd(q_at_lens,el.abcd);
+            x_last_calc = el.x;
+        end
+    end
+    
+    %points to the right of the last optical element
+    index = (x_last_calc <= x_pos);
+    x=x_pos(index);
+    q(index)= q_after_free_space(q_last_calc,x-x_last_calc);
+       
+	
+	
 end
 
 %!test