%--------------------------------------------------------------------------
% [Temp] =
% FT_temperature_field_from_coating_absorption(beam,mirror,r,z,num_terms)
%
% A function for Matlab which approximates the temperature field in a
% mirror heated by a beam as a sum of bessel functions, where the heating
% is due to absorption in the mirror coating.
%
% N.B. currently this function is only valid for an LG00 incident 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]
% z:        Depth coordinate along beam axis [m]
% num_terms:Number of terms to use in approximation
%
% Temp:     Returned temperature field
%
% Equations taken from Vinet 2009, Thermal issues in Advanced GW detectors 
% and Virgo Optics book.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    07.02.2012
%--------------------------------------------------------------------------
%

function [Temp] = FT_temperature_field_from_coating_absorption(beam,mirror,r,z,num_terms)

    % Currently only valid for LG00
    p = 0;
    l = 0;
    
    % Beam and mirror parameters
    P = beam.P;
    w = beam.w;
    Cabs = mirror.Cabs;
    a = mirror.a;
    h = mirror.h;
    K = mirror.K;
    X = mirror.X;
    
    % Beam is axis symmetric
    r_abs = abs(r);
    
    % 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);
    gamma_0s = zeta_0s*h/(2*a);
    d_10s = zeta_0s.*sinh(gamma_0s) + X*cosh(gamma_0s);
    d_20s = zeta_0s.*cosh(gamma_0s) + X*sinh(gamma_0s);

    % Approximating temperature field as bessel sum
    for j=1:length(z)
        T = 0;
        for s=1:num_terms
            T_term = Cabs*P/(2*pi*K*a) * p_0s(s) * ( cosh(zeta_0s(s)*z(j)/a)/d_10s(s) - sinh(zeta_0s(s)*z(j)/a)/d_20s(s));  
            T = T + T_term*besselj(0,zeta_0s(s)*r_abs/a);
        end
        Temp(j,:) = T(:);
    end

end

