function Gouy_phase = Calculate_Gouy_Phase(Ein,z)

[family,m,n] = Read_mode_name(Ein.Mode_name);
HOM_order = append(family,' ' ,num2str(n+1),' ' ,num2str(m));

if n==0 && m==0

    phi_PW = 2*pi/Ein.Wavelength *z;

    Ein_z = Propagate_E(Ein,z);
    E_overlap = Ein.Field.*conj(Ein_z.Field); 

    Gouy_phase_matrix = atan2(imag(E_overlap),real(E_overlap))-phi_PW;
    Gouy_phase = rem(Gouy_phase_matrix(Ein.Grid.Num_point/2,Ein.Grid.Num_point/2),2*pi);

%     figure(1)
%     E_Plot(Ein_z); drawnow

else

%     TO BE DONE

end