Gaussian mode

Return a Gaussian mode as a 1xP vector l: mode number+1 p: number of modes R: end radius

function mode = GaussianMode(l,r,z,lambda,omega0,n)

    %Constants

    k       = 2*pi * n/ lambda;
    z0      = pi * omega0^2 * n / lambda;
    omega   = omega0 * sqrt(1 + (z /z0).^2 );

    %Find the mode vector
    if(z==0)
        expled  = 1;
    else
        R       = z * (1 + (z0 / z).^2 );
        expled  = exp(1i * k * r.^2 / ( 2 * R ) );
    end
    norm    = sqrt(2/pi) / omega0;
    mode    = norm .* omega0 ./ omega .* exp( - r.^2 ./ omega^2 ) .* polyval(LaguerreGen(l-1, 0), 2*r.^2 ./ (omega.^2 ) ) .* exp( (2*(l-1)+1)*1i * atan( z / z0) ) .* expled ;

end