%--------------------------------------------------------------------------
% function [Z_coat,coeffs] = FT_thermal_lens_from_coating_absorption(beam,mirror,r,num_terms)
%
% A function for Matlab which approximates the thermal lens resulting from 
% a beam incidnet on a mirror described by a sum of bessel functions.  The 
% steady state solution is taken and the thermal lens is calculated due to 
% absorption in the mirror coating.
%
% N.B. Currently only valid for an incidnet LG00 beam.
% 
% beam:     Structure storing beam parameters (see FT_init_thermal_beam.m)
% mirror:   Structure storing mirror parameters (see
%           FT_init_thermal_mirror.m)
% r:        Radial coordinate [m]
% num_terms:Number of terms to use in approximation. (50 is usually
%           sufficient).
%
% Z_coat:   Approximate thermal lens [m]
% coeffs:   Coefficients used to approximate thermal lens from sum of
%           bessel functions.
%
% Equations taken from Vinet 2009, Thermal issues in GW detectors and the
% Virgo Optics book.
% 
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    09.02.2012
%--------------------------------------------------------------------------
%

function [Z_coat,coeffs] = FT_thermal_lens_from_coating_absorption(beam,mirror,r,num_terms)
    
    % Mirror and beam parameters 
    % Currently only valid LG00
    p = 0;
    l = 0; 
    a = mirror.a;
    h = mirror.h;
    K = mirror.K;
    dn_dT = mirror.dn_dT;
    Cabs = mirror.Cabs;
    X = mirror.X
    P = beam.P;
    lambda = beam.lambda;
    w = beam.w;
    
    % Lens is axis-symmetric
    r_abs = abs(r);

    % Constant
    c = dn_dT * Cabs*P / (pi*K); 
    
    % Calculating coefficients
    zeta_0s = FT_bessel_zeros(0,X,num_terms);
    p_0s = FT_calculate_fourier_bessel_coefficients(0,0,X,w,a,num_terms);
    
    % Calculating sum of terms to approximate thermal lens
    Z_coat=0;
    
    for s=1:num_terms
        
        coeffs(s) = c*p_0s(s)*(1/zeta_0s(s)) * (1-exp(-zeta_0s(s)*h/a))/(zeta_0s(s)+X-(zeta_0s(s)-X)*exp(-zeta_0s(s)*h/a));
        Z_coat = Z_coat +  coeffs(s)*besselj(0,zeta_0s(s)*r_abs/a);
    end
          

end

