%--------------------------------------------------------------------------
% function [uz_coat,duz_coat] =
% FT_thermal_distortion_from_coating_absorption(beam,mirror,r,num_terms)
%
% A function for Matlab which calculates an approximation of the distortion
% of a mirror due to heating by an incident laser beam and absorption in 
% the mirror coating, in the steady state.
%
% N.B. Currently only valid for an incidnet LG00 beam.
%
% beam:     Structure storing the parameters of the beam (see
%           FT_init_thermal_beam.m).
% mirror:   Structure storing thermal properties for the mirror (see
%           FT_init_thermal_mirror.m).
% r:        Radial coordinate [m]
% num_terms:Number of terms to include in approximation.  (50 terms is
%           usually sufficient).
% 
% U_coat:  Returned distortion of the mirror surface [m]
%
% 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    11.02.2012
%--------------------------------------------------------------------------
%

function [U_coat] = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,num_terms)

    % Mirror and beam parameters 
    % Currently only valid for LG00
    p = 0;
    l = 0;  
    
    % Mirror parameters
    a = mirror.a;
    h = mirror.h;
    K = mirror.K;
    Cabs = mirror.Cabs;
    X = mirror.X;
    alpha = mirror.alpha;
    Y = mirror.Y;
    pois = mirror.pois;
    
    % Beam parameters
    P = beam.P;
    lambda = beam.lambda;
    w = beam.w;
    
    % Lens is axis-symmetric
    r_abs = abs(r);
    
    % Constants
    c1 = alpha*(1+pois)*Cabs*P/(pi*K);
    c2 = -12*alpha*Y*Cabs*P*X/(pi*a*K*h^3);
    
    % Calculate coefficients
    zeta_0s = FT_bessel_zeros(0,X,num_terms);
    
    % Approximation of surface distortion
    uz_coat = 0;
    B = 0;
    
    for s=1:num_terms
        
        f1 = c1 * zeta_0s(s)*exp(-zeta_0s(s)^2*w^2/(8*a^2)) / ((zeta_0s(s)^2+X^2)*besselj(0,zeta_0s(s))^2) * (zeta_0s(s)+X-(zeta_0s(s)-X)*exp(-2*zeta_0s(s)*h/a)) / ((zeta_0s(s)+X)^2-(zeta_0s(s)-X)^2*exp(-2*zeta_0s(s)*h/a)); 
        uz_coat = uz_coat + f1 * (besselj(0,zeta_0s(s)*r_abs/a)-1);
        
        % Coefficient for Saint-Venant correction
        f2 = c2 * exp(-w^2*zeta_0s(s)^2/(8*a^2)) / ((zeta_0s(s)^2+X^2)*besselj(0,zeta_0s(s))) * (1/((zeta_0s(s)+X)^2-((zeta_0s(s)-X)^2*exp(-2*zeta_0s(s)*h/a))));
        f2 = f2 * ( ((zeta_0s(s)-X)*exp(-zeta_0s(s)*h/a)*a^2/(zeta_0s(s)^2)* ( (zeta_0s(s)*h/(2*a) * (exp(zeta_0s(s)*h/(2*a)) + exp(-zeta_0s(s)*h/(2*a)))) + (exp(-zeta_0s(s)*h/(2*a))-exp(zeta_0s(s)*h/(2*a)))) + ((zeta_0s(s)+X)*a^2/(zeta_0s(s)^2)*(-(zeta_0s(s)*h/(2*a)*(exp(-zeta_0s(s)*h/(2*a))+exp(zeta_0s(s)*h/(2*a)))) + (exp(zeta_0s(s)*h/(2*a))-exp(-zeta_0s(s)*h/(2*a)))))));
                
        B = B + f2;
       
    end
    
    % Saint-Venant correction
    duz_coat = (1-pois)/(2*Y) * B * r.^2;
    
    % Final thermal distortion
    U_coat = uz_coat + duz_coat;

end

