trial_solution;

q0=wr2q(w0,r0,lambda);
x=linspace(x0,xf,1000); % we will calculate beam profile between x0 and xf
fprintf('======== Forward propagation ======\n')
q_forward=gbeam_propagation(x,q0,x0,optics);
[w_forward,r_forward]=q2wr(q_forward, lambda);

fprintf('======== Backward propagation =====\n')
q_backward0=wr2q(wf,rf,lambda);

q_backward=gbeam_propagation(x,q_backward0,xf,optics);
[w_backward,r_backward]=q2wr(q_backward, lambda);

fprintf('========= Plotting ================\n')
plot ( ...
	x,w_forward, '-r',  ...
	x,-w_forward, '-r', ...
	x, w_backward, '-.b', ...
	x, -w_backward, '-.b')
legend({'forward propagation', '', 'backward propagation', ''})      

fprintf('======= final check =============\n')
fprintf('following are desired values: \n')
w0
r0
x0
wf
rf
xf
fprintf('======= after propagation ========\n')
fprintf('values below should match wf and rf: \n')
[waste_at_the_end,radius_at_the_end]  = q2wr(q_forward(end), lambda);
waste_at_the_end
radius_at_the_end