function fs_abcd = abcd_free_space( distance)
	fs_abcd=[1, distance; 0,1];
endfunction

function lens_abcd =abcd_lens(focal_distance)
	lens_abcd = [1, 0; -1/focal_distance, 1];
endfunction

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 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_forward(x_pos, q_in, x_in,  optics_elements)
% calculate the 'q' parameter of the Gaussian beam propagating through optical
% 'optics_elements' only in the positive direction  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
% x_pos must be monotonic!
	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(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_afteer_element(q_last_calc,abcd);
				q_last_calc=q_afteer_element(q_last_calc,el.abcd);
				x_last_calc=el.x;
			endif
		endfor
		if (x_pos_i > x_last_calc);
			abcd=abcd_free_space(x_pos_i-x_last_calc);
			q_last_calc=q_afteer_element(q_last_calc,abcd);
			x_last_calc=x_pos_i;
		endif
		q(i)=q_last_calc;
	endfor
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
% x_pos must be monotonic!

	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)
	for i=1:size(q,2)
		waste(i)=sqrt (-lambda/pi/imag(1/q(i)));
	endfor
endfunction

function radius =q2radius(q, lambda)
	for i=1:size(q,2)
		radius(i)=(1/real(1/q(i)));
	endfor
endfunction

function q=waste_r2q(waste,R,lambda)
	q=1/(1/R-I*lambda/pi/waste/waste);
endfunction