function [q] = gbeam_propagation_froward_only(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
% Note!: if optics element position coinside with a x_pos at some place
%   the reported value will be the one calculated after the element
%   so the resulting value will coinside with submitted q_in
%
% all x_pos must be to the right of x_in
% x_pos must be monotonic!
    
    %Initialize q_lens (q at position of lens)
    q_lens = zeros(1,size(optics_elements,2));

	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
	q=0*x_pos; % q vector initialization
	q_last_calc=q_in;
	x_last_calc=x_in; % the furthest calculated point
    
    
    for k = 1:length(optics_elements)
        el=optics_elements{k};
        elx = el.x;
        if elx > x_last_calc
            %Find all points to the left of the current optical element
            index = (x_last_calc <= x_pos) & (x_pos < elx);
            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,elx - x_last_calc);
            %Applying lens transformation
            q_last_calc=q_after_abcd(q_at_lens,el.abcd);
            x_last_calc = elx;
        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 
%! 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.8098e-96 + 1.3453e+02i;
%! x0=0;
%! optics={lns1,lns2,lns3};
%! x=[0, 0.1, 0.2, 0.3, 0,4,  0.5, 0.6, 0.7];
%! q_test = [1.8098e-96 + 1.3453e+02i,  1.0000e-01 + 1.3453e+02i,  -1.5260e-01 + 1.7310e-04i, -5.2600e-02 + 1.7310e-04i,  -5.2600e-02 + 1.7310e-04i,  3.4371e+00 + 1.8961e-03i, 3.4371e+00 + 1.8961e-03i,  3.4371e+00 + 1.8961e-03i,  3.4371e+00 + 1.8961e-03i];
%! 
%! q=gbeam_propagation(x,q0,x0,optics);
%! assert(q,q_test, -1e-4)