function [q, q_lens] = gbeam_propagation(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), q_lens] = gbeam_propagation_froward_only(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=fliplr(optics_elements);
		% we need to flip all optics elements around x_in as well
        n_elements = length(optics_elements);
        
		for i=1:n_elements
			optics_elements_backw{i}.x=x_in-optics_elements_backw{i}.x;
            % there is a tricky case: if position of any optics element coincides
			% with x_in i.e it is 0 in backwards x coordinates (or at index 1) then
			% we need to apply abcd matrix of that element in advance.
			if ( optics_elements_backw{i}.x == 0) 
				q_in_backw(1)  = q_after_abcd( q_in_backw(1), optics_elements_backw{i}.abcd);
			end

		end


		q_backw = gbeam_propagation_froward_only(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

    
    
end

%!test 
%! lambda= 1.064E-6 ;
%! f1=0.1526;
%! lns1.abcd=abcd_lens( f1) ;
%! lns1.x= 0.2;
%! f2=0.019;
%! lns2.abcd=abcd_lens( f2) ;
%! lns2.x= lns1.x+ f1 + f2;
%! lns3.abcd=abcd_lens( 0.0629) ;
%! lns3.x= 0.500;
%! q0=-1.5260e-01 + 1.7310e-04i;
%! x0=0.2;
%! optics={lns1,lns2,lns3};
%! x=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7];
%! q0=1.8098e-96 + 1.3453e+02i;
%! x0=0;
%! q_test = gbeam_propagation(x,q0,x0,optics);
%! 
%! % test beam is the same but starting point is taken midway 
%! % and coincides with optics element{1}.x since (x(3) = 0.2)
%! q=gbeam_propagation(x,q_test(3),x(3),optics);
%! assert(q,q_test, -1e-6)