%--------------------------------------------------------------------------
% FTE_aLIGO_thermal_lenses.m
%
% Matlab example which calculates the thermal distortions on the input test 
% mass (itm) and end test mass (etm) in an Advanced LIGO arm cavity with 
% 1ppm absorption in the mirrors.  This example also produces finesse map 
% files of the thermal distortions for use in simulations.  It also returns
% the radii of curvature for the distortions.
%
% In this example the overall thermal distortion is calculated and 
% approximated by a parabalid (i.e. curved surface).  The residual
% distortion (overal - approximate parabaloid) is also calculated.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    28.01.2013
%--------------------------------------------------------------------------
%

clear all;

% Choose coating absorption
option = 'coating';
    
% Filenames for map files
file1 = 'etm_thermal_lens.txt';
file2 = 'etm_thermal_res_lens.txt';
file3 = 'itm_thermal_lens.txt';
file4 = 'itm_thermal_res_lens.txt';
filea1 = 'etm_aperture.txt';
filea2 = 'itm_aperture.txt';

% Number
num_terms = 100;

% Beam parameters (aLIGO parameters)
% Two beam parameters at the two cavity mirrors
P = 770e3;
p = 0;
l = 0;
w_itm = 5.3e-2;
w_etm = 6.2e-2;
lambda = 1064e-9;
    
% Initialising thermal beam structures
beam_etm = FT_init_thermal_beam(P,p,l,w_etm,lambda);
beam_itm = FT_init_thermal_beam(P,p,l,w_itm,lambda);
    
% Mirror material parameters
% Case for 1ppm absorption
Cabs = 1e-6;
K = 1.38;
E = 1;
alpha = 5.4e-7;
Y = 7.3e10;
pois = 0.17;
a = 0.17;
h = 0.2;
dn_dT = -1.1e-5;
% For 1 ppm absorption
Babs = 1e-6/(P*h);
    
% Intiate thermal mirror
mirror = FT_init_thermal_mirror(dn_dT,Cabs,Babs,K,E,alpha,Y,pois,a,h);
    
% Radial coordinate
r = -a:2*a/256:a;
    
% Thermal lens: etm
[Z_etm,c_etm,p_etm] = FT_thermal_lens(beam_etm,mirror,r,option,num_terms);

% Lens described by approximate parabaloid: etm 
Zp_etm = c_etm*r.^2 + p_etm;

% Residual lens: etm
Zr_etm = Z_etm - Zp_etm;

% Creating mirror maps representing overall and residual lenses: etm    
[Zmap_etm,rmap_etm,Rc_etm] = FT_phase_map_from_thermal_lens(beam_etm,mirror,r,option,num_terms);
 
% Creating aperture map and writing maps to file: etm
Amap_etm = FT_create_aperture_map(Zmap_etm);
FT_write_surface_map(file1,Zmap_etm);
FT_write_surface_map(file2,rmap_etm);
FT_write_surface_map(filea1,Amap_etm);
    
% Thermal lens: itm
[Z_itm,c_itm,p_itm] = FT_thermal_lens(beam_itm,mirror,r,option,num_terms);

% Lens described by approximate parabaloid: itm
Zp_itm = c_itm*r.^2 + p_itm;

% Residual lens: itm
Zr_itm = Z_itm - Zp_itm;
    
% Creating mirror maps representing overall and residual lenses: itm
[Zmap_itm,rmap_itm,Rc_itm] = FT_phase_map_from_thermal_lens(beam_itm,mirror,r,option,num_terms);
    
% Creating aperture maps and writing maps to file: itm
Amap_itm = FT_create_aperture_map(Zmap_itm);
FT_write_surface_map(file3,Zmap_itm);
FT_write_surface_map(file4,rmap_itm);
FT_write_surface_map(filea2,Amap_itm);
    
% Plotting results
figure
    % Plotting etm thermal lens and approximate parabaloid
    hold on;
    plot(r,Z_etm,'r','LineWidth',3)
    plot(r,Zp_etm,'b--','LineWidth',3)
    grid on;
    xlabel('r [m] ')
    ylabel('Thermal lens [m] ')
    title('Thermal distortion: ETM')
    legend('Thermal lens','Approximate paraboloid')
    set(gca,'XLim',[-a a])
                
figure
    % Plotting etm residual lens
    plot(r,Zr_etm,'r','LineWidth',3)
    grid on;
    xlabel('r [m] ')
    ylabel('Thermal lens [m] ')
    title('Residual thermal lens: ETM')
    set(gca,'XLim',[-a a])

% Plotting thermal distortion as map: etm
FT_plot_mirror_map(Zmap_etm,0)

% Plotting residual distortion as map: etm
FT_plot_mirror_map(rmap_etm,0)
        
figure
    % Plotting itm thermal distortion and approximate parabaloid
    hold on;
    plot(r,Z_itm,'r','LineWidth',3)
    plot(r,Zp_itm,'b--','LineWidth',3)
    grid on;
    xlabel('r [m] ')
    ylabel('Thermal lens [m] ')
    title('Thermal lens: ITM')
    legend('Thermal lens','Approximate paraboloid')
    set(gca,'XLim',[-a a])
                
figure
    % Plotting itm residual distortion
    plot(r,Zr_itm,'r','LineWidth',3)
    grid on;
    xlabel('r [m] ')
    ylabel('Thermal lens [m] ')
    title('Residual thermal lens: ITM')
    set(gca,'XLim',[-a a])

% Plotting thermal lens as map: itm
FT_plot_mirror_map(Zmap_itm,0)
    
% Plotting residual lens as map: itm
FT_plot_mirror_map(rmap_itm,0)
    
% Outputting radii of curvature of lenses (of approximate parabaloid)
disp(sprintf('ETM approximate (gaussian weighted) radius of curvature: %g [m]',Rc_etm))
disp(sprintf('ITM approximate (gaussian weighted) radius of curvature: %g [m]',Rc_itm))



