lot's of changes to make forward and backward propagation user independent

Ignore-this: d9d6aa99367b8f2f682338722aea4c9c

darcs-hash:20110413054232-067c0-dd341b328933a7e7e06375ff5339077e4f585e32

1 files changed, 75 insertions, 10 deletions

diff --git a/abcd.m b/abcd.m
index 8404314..7e02200 100644
--- a/abcd.m
+++ b/abcd.m
@@ -10,23 +10,30 @@ function qnew=q_afteer_element(q_old,abcd)
 	qnew=(q_old*abcd(1,1)+abcd(1,2))/(q_old*abcd(2,1)+abcd(2,2));
 endfunction
 
-function q = prop(x_pos, q_in, elements)
+function q = prop_forward(x_pos, q_in, x_in,  optics_elements)
 % calculate the 'q' parameter of the Gaussian beam propagating through optical
-% 'elements' array along 'x' axis at points 'x_pos'
-%  takes the gaussian beam with initial q_in parameter at x_pos(1)
+% 'optics_elements' array along 'x' axis at points 'x_pos'
+%  takes the gaussian beam with initial q_in parameter at x_in
+%
+% all x_pos must be to the right of x_in
+	if (any(x_pos < x_in))
+		error('all beam positions must be to the right of the x_in');
+	end
 
+	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(elements); 
+	Nel=length(optics_elements) ;
 	q=0*x_pos; % q vector initialization
-	q(1)=q_in;
 	q_last_calc=q_in;
-	x_last_calc=x_pos(1); % the furthest calculated point
-	for i=2:Np
+	x_last_calc=x_in; % the furthest calculated point
+	for i=1:Np
 		x_pos_i=x_pos(i);
-		for k=1:length(elements) 
-			% iterates through optical elements to make sure 
+		for k=1:length(optics_elements) 
+			% iterates through optics_elements to make sure 
 			% we take them in account for the beam propagation
-			el=elements{k};
+			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_afteer_element(q_last_calc,abcd);
@@ -41,6 +48,64 @@ function q = prop(x_pos, q_in, elements)
 		endif
 		q(i)=q_last_calc;
 	endfor
+end
+
+function optics = arrange_optics_along_x(optics_unsorted)
+% arrange optics in proper order so it x position increases with number
+	N=length(optics_unsorted);
+
+	% assign x positions
+	x=zeros(1,N);
+	for i=1:N
+		x(i)=optics_unsorted{i}.x;
+	end
+
+	[xs,indx]=sort(x);
+	cntr=1;
+	for i=indx
+		optics{cntr}=optics_unsorted{i};
+		cntr=cntr+1;
+	end
+end
+
+function q = prop(x_pos, q_in, x_in,  optics_elements)
+% calculate the 'q' parameter of the Gaussian beam propagating through optical
+% 'optics_elements' array along 'x' axis at points 'x_pos'
+%  takes the gaussian beam with initial q_in parameter at x_in
+
+	q=0*x_pos; % q vector initialization
+
+	if any(x_pos >= x_in)
+		% Forward propagation to the right of x_in
+		q(x_pos >= x_in) = prop_forward(x_pos(x_pos>=x_in), q_in, x_in,  optics_elements);
+	end
+
+	if any(x_pos < x_in)
+		% Backward propagation part the left of x_in
+		% do it as forward propagation of the reverse beam
+		x_backw=x_pos(x_pos<x_in);
+		% now let's reflect the beam with respect to x_in
+		% and solve the problem as forward propagating.
+		x_backw=x_in-x_backw;
+		% now we need to flip x positions
+		x_backw=fliplr(x_backw);
+		% reflected beam means inverted radius of curvature or real part of q parameter
+		q_in_backw = -real(q_in) + 1i*imag(q_in);
+		optics_elements_backw=optics_elements;
+		% we need to flip all optics elements around x_in as well
+		for i=1:length(optics_elements_backw)
+			optics_elements_backw{i}.x=x_in-optics_elements_backw{i}.x;
+		end
+
+		q_backw = prop_forward(x_backw, q_in_backw, 0,  optics_elements_backw);
+		% now we need to flip the radius of curvature again
+		q_backw =  -real(q_backw) + 1i*imag(q_backw);
+
+		% final assignment of the backwards propagating beam 
+		% which we need to flip back
+		q(x_pos<x_in) = fliplr(q_backw);
+	end
+
 endfunction
 
 function waste =q2waste(q, lambda)