diff --git a/src/new/thermal/FT_thermal_distortion.m b/src/new/thermal/FT_thermal_distortion.m
index 3ce9b81bf95ed520cdd54f0e21393fef7f156d60..f98b18e008eeadd32fa233b50c57b194b6a8c35c 100644
--- a/src/new/thermal/FT_thermal_distortion.m
+++ b/src/new/thermal/FT_thermal_distortion.m
@@ -6,8 +6,6 @@
% The function returns the distortion of the surface and the parameters c
% and p for the closest (Gaussian weighted) fit paraboloid.
%
-% N.B. Currently only valid for coating absorption.
-%
% beam: Structure storing beam parameters (FT_init_thermal_beam.m)
% mirror: Structure storing mirror parameters (FT_init_thermal_mirror.m)
% r: Radial coordinate [m]
@@ -29,13 +27,14 @@ function [U,c,p] = FT_thermal_distortion(beam,mirror,r,option,num_terms)
% Calculate thermal distortion from coating absorption
if strcmp(option,'coating')
- U = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,num_terms);
- [c,p] = FT_approximate_paraboloid_for_thermal_distortion(beam,mirror,num_terms);
+ [U,coeffs,SV] = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,num_terms);
+ [c,p] = FT_approximate_paraboloid_for_thermal_distortion(coeffs,SV,beam,mirror);
% Calculate distortion from bulk absorption
elseif strcmp(option,'bulk')
- % Currently only valid for coating absorption
+ [U,coeffs] = FT_thermal_distortion_from_bulk_absorption(beam,mirror,r,num_terms);
+ [c,p] = FT_approximate_paraboloid_for_thermal_distortion(coeffs,0,beam,mirror);
end
diff --git a/src/new/thermal/FT_thermal_distortion_from_bulk_absorption.m b/src/new/thermal/FT_thermal_distortion_from_bulk_absorption.m
new file mode 100644
index 0000000000000000000000000000000000000000..6dd46b0e6bab367a29618326d3bf725d676dc2ed
--- /dev/null
+++ b/src/new/thermal/FT_thermal_distortion_from_bulk_absorption.m
@@ -0,0 +1,80 @@
+%--------------------------------------------------------------------------
+% function [U_bulk,coeffs] = FT_thermal_distortion_from_bulk_absorption(beam,mirror,r,num_terms)
+%
+% A Matlab function which calculates an approximation of the thermal
+% distortion of a mirror surface due to absorption of power in the bulk of
+% the mirror (or substrate). The approximation expands the distortion as a
+% series of bessel functions, for a mirror in equilibrium. The equations
+% are taken from the Living Review, 'On Special Optical Modes and Thermal
+% Issues in Advanced GW Interferometric Detectors', Vinet, 2009.
+%
+% N.B. Currently only valid for an incident LG00 beam.
+%
+% beam: Structure containing information about the incident beam (see
+% FT_init_thermal_beam.m).
+% mirror: Structure containing thermal properties of the mirror (see
+% FT_init_thermal_mirror.m).
+% r: Vector of radial coordinates [m]
+% num_terms:Number of bessel terms to include in series.
+%
+% U_bulk: Returned surface distortion [m]
+% coeffs: Coefficients for bessel series expansion.
+%
+% Part of the Simtools package, http://www.gwoptics.org/simtools
+% Charlotte Bond 13.11.2013
+%--------------------------------------------------------------------------
+%
+
+function [U_bulk,coeffs] = FT_thermal_distortion_from_bulk_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;
+ Babs = mirror.Babs;
+ X = mirror.X;
+ alpha = mirror.alpha;
+ Y = mirror.Y;
+ pois = mirror.pois;
+
+ % Beam parameters
+ Power = beam.P;
+ lambda = beam.lambda;
+ w = beam.w;
+
+ % Lens is axis-symmetric
+ r_abs = abs(r);
+
+ % Constants
+ C1 = alpha*(1+pois)*Babs*Power*a/(pi*K);
+ zeta_n = FT_bessel_zeros(l,X,num_terms);
+ pn = FT_calculate_fourier_bessel_coefficients(p,l,X,w,a,num_terms);
+ gamma_n = zeta_n*h/(2*a);
+ Tn = sinh(gamma_n).*cosh(gamma_n)+gamma_n;
+ d1n = zeta_n.*sinh(gamma_n)+X*cosh(gamma_n);
+
+ U_bulk = 0;
+
+ % Calculate coefficents for thermal distortion in terms of bessel
+ % functins
+ for n=1:num_terms
+
+ C2 = pn(n)/zeta_n(n)^3 * sinh(gamma_n(n));
+ t1 = 2*sinh(gamma_n(n))/Tn(n);
+ t2 = X/d1n(n);
+
+ J0 = besselj(0,zeta_n(n)*r_abs/a);
+
+ coeffs(n) = -C1*C2*(t1-t2);
+
+ U_bulk = U_bulk + coeffs(n)*(1-J0);
+
+ end
+
+end
+
diff --git a/src/new/thermal/FT_thermal_distortion_from_coating_absorption.m b/src/new/thermal/FT_thermal_distortion_from_coating_absorption.m
index 4dbb2d538357282b6bcab0fcaaaf47e5acf4fff7..1f4a6112cb492c1585e36fedffa9a49edb275603 100644
--- a/src/new/thermal/FT_thermal_distortion_from_coating_absorption.m
+++ b/src/new/thermal/FT_thermal_distortion_from_coating_absorption.m
@@ -16,7 +16,10 @@
% num_terms:Number of terms to include in approximation. (50 terms is
% usually sufficient).
%
-% U_coat: Returned distortion of the mirror surface [m]
+% U_coat: Returned distortion of the mirror surface [m]
+% coeffs: Coefficients in expansion of thermal distortion as sum of bessel
+% functions.
+% SV: Coefficients for Saint-Venant correction.
%
% Equations taken from Vinet 2009, Thermal issues in GW detectors and the
% Virgo Optics book.
@@ -26,7 +29,7 @@
%--------------------------------------------------------------------------
%
-function [U_coat] = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,num_terms)
+function [U_coat,coeffs,SV] = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,num_terms)
% Mirror and beam parameters
% Currently only valid for LG00
@@ -44,7 +47,7 @@ function [U_coat] = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,
pois = mirror.pois;
% Beam parameters
- P = beam.P;
+ Power = beam.P;
lambda = beam.lambda;
w = beam.w;
@@ -52,34 +55,36 @@ function [U_coat] = FT_thermal_distortion_from_coating_absorption(beam,mirror,r,
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);
+ c1 = alpha*(1+pois)*Cabs*Power/(2*pi*K);
+ c2 = -12*alpha*Y*X*Cabs*a*Power/(pi*K*h^3);
% Calculate coefficients
- zeta_0s = FT_bessel_zeros(0,X,num_terms);
+ zeta_n = FT_bessel_zeros(2*l,X,num_terms);
+ pn = FT_calculate_fourier_bessel_coefficients(p,l,X,w,a,num_terms);
+ gamma_n = zeta_n*h/(2*a);
+ d1n = zeta_n.*sinh(gamma_n)+X*cosh(gamma_n);
+ d2n = zeta_n.*cosh(gamma_n)+X*sinh(gamma_n);
% Approximation of surface distortion
- uz_coat = 0;
- B = 0;
+ uz = 0;
+ w1 = 0;
+ % Calculate distortion as bessel sum
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;
+ % Negative compared to Vinet 2009, for match Finesse convention
+ coeffs(s) = - c1 * pn(s)/zeta_n(s) * (sinh(gamma_n(s))/d1n(s) + cosh(gamma_n(s))/d2n(s));
+ w1 = w1 - c2*pn(s)*besselj(0,zeta_n(s))/(zeta_n(s)^4)*(gamma_n(s)*cosh(gamma_n(s))-sinh(gamma_n(s)))/d2n(s);
+ uz = uz + coeffs(s)*(1-besselj(2*l,zeta_n(s)*r_abs/a));
end
% Saint-Venant correction
- duz_coat = (1-pois)/(2*Y) * B * r.^2;
+ SV = -(1-pois)/(2*Y)*w1;
+ duz = SV * r.^2;
- % Final thermal distortion
- U_coat = uz_coat + duz_coat;
+ % Final thermal distortion
+ U_coat = uz + duz;
end