%--------------------------------------------------------------------------
% FTE_virgo_thermal_lens.m
%
% A Matlab example which calulates the thermal lens and temperature field 
% in a mirror in the arm cavities of the virgo detector.  This example 
% plots the thermal lens, it's approximate parabaloid, the residual
% (difference from the approximate parabaloid) and the temperature field in
% the mirror.  You can compare the results with pages 289 - 301 in the 
% Virgo Optics Book.  The thermal lens is calculated from 1W absorption in 
% the coating.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    09.02.2012
%--------------------------------------------------------------------------
%

clear all;

% Coating absorption
option = 'coating';

% Number of terms to include in the expansion
num_terms = 100;

% Beam parameters (Virgo parameters)
P = 1;
p = 0;
l = 0;
w = 2e-2;
lambda = 1064e-9;
    
% Initiate thermal beam
beam = FT_init_thermal_beam(P,p,l,w,lambda);
    
% Mirror parameters
dn_dT = -0.87e-5;
% For 1W absorption have Cabs = 1 and P = 1
Cabs = 1;
K = 1.38;
E = 0.8;
a = 0.175;
h = 0.1;
% For 1W absorption
Babs = 1/(P*h);
    
% Initiate thermal mirror
mirror = FT_init_thermal_mirror(dn_dT,Cabs,Babs,K,E,0,0,0,a,h);
    
% Radial and z (depth) coordinates
r = -0.06:0.001:0.06;    
z = -0.05:0.001:0.05;
    
% Thermal lens
[z_lens,c,p] = FT_thermal_lens(beam,mirror,r,option,num_terms);

% Approximate parabaloid
z_para = c*r.^2 + p;

% Residual thermal lens (difference from approximate parabaloid)
z_res = z_lens - z_para;
    
% Temperature field
T = FT_temperature_field(beam,mirror,r,z,option,num_terms);
    
% Plot results
% Plotting the thermal lens and approximate parabaloid
figure 
    hold on;
    plot(r,z_lens,'r','LineWidth',3)
    plot(r,z_para,'b--','LineWidth',3)
    grid on;
    xlabel('r [m] ')
    ylabel('Thermal lensing [m/W] ')
    title('Thermal lensing')
    set(gca,'XLim',[-0.06 0.06])
    legend('Thermal lens','Approximate paraboloid')

% Plotting the residual thermal lens
figure 
    hold on;
    plot(r,z_res,'r','LineWidth',3)
    grid on;
    xlabel('r [m] ')
    ylabel('Thermal lens [m/W] ')
    set(gca,'XLim',[-0.04 0.04])
    title('Thermal lens: difference from approximae parabaloid ')
      
% Plotting the temperature field
figure
    surfc(r,z,T,'edgecolor','none')
    axis tight;
    %set(gca,'XLim',[-0.05 0.05])
    %set(gca,'YLim',[-0.05 0.05])
    title('Temperature field ')
    xlabel('r [m] ')
    ylabel('z [m] ')
    hc=colorbar;
    ylabel(hc,'Temperature [K/W] ');


