%--------------------------------------------------------------------------
% function [I] = FT_FB_expansion_of_LG_beam(p,l,w,a,r,X,num_terms)
%
% A function for Matlab which calculates the intensity distribution of a
% Laguerre-Gauss beam in terms of an expansion in a Fourier-Bessel (FB)
% series.  The intesnity is the normalised intensity, i.e. 1W of power in
% the beam.
%
% p,l:      Indices of the (helical) LG mode
% w:        Beam radius [m]
% a:        Mirror radius [m]
% r:        Radial coordinates [m]
% X:        Reduced radiation constant (X = 4*o_*T0^3*a/K) where o_ is SB
%           constant corrected for emmistivity (o_ typically 0.8-0.9 o_SB)
%           T0 = temperature of vaccum walls, K = thermal conductivity of
%           mirror.
% num_terms:The number of terms to include in the expansion.  (N.B. A good
%           recreation is generally achieved with 50 terms).
%
% I:        Returned intensity.
%
% Equations taken from Vinet 2009, Thermal Issues in Advanced GW detectors 
% and Virgo book on optics.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond 07.02.2012
%--------------------------------------------------------------------------
%

function [I] = FT_FB_expansion_of_LG_beam(p,l,w,a,r,X,num_terms)
    
    % I is the normalised beam.
    P=1;    
    
    % Solutions of bessel equation y*J'_n(y)+X*J_n(y)=0
    zeta_ls = FT_bessel_zeros(2*l,X,num_terms);
    
    % Calculate fourier bessel coefficients to expand the beam in terms of
    % bessel functions
    p_ls = FT_calculate_fourier_bessel_coefficients(p,l,X,w,a,num_terms);
    
    % Intensity is axis-symmetric
    r_abs = abs(r);
    
    I = 0;
    
    % Compute I as sum of bessel terms
    for s=1:num_terms
        
        I = I + P/(pi*a^2)*p_ls(s)*besselj(2*l,zeta_ls(s)*r_abs/a);        
        
    end

    
end

