diff --git a/src/new/thermal/FT_thermal_lens.m b/src/new/thermal/FT_thermal_lens.m
index 64d0e6e6ac338a96051af3c6c70685a7dabcb20f..d9d4a51430bf9e80674f03905355d04bd1b3d162 100644
--- a/src/new/thermal/FT_thermal_lens.m
+++ b/src/new/thermal/FT_thermal_lens.m
@@ -6,6 +6,8 @@
 % returns the lens as well as the approximate (Gaussian weighted) 
 % paraboloid.
 %
+% N.B. Currently only valid for an incident LG00 mode.
+%
 % beam:     Structure containing beam parameters (FT_init_thermal_beam.m)
 % mirror:   Structure containing mirror parameters 
 %           (FT_init_thermal_mirror.m)
@@ -32,7 +34,7 @@ function [Z,c,p] = FT_thermal_lens(beam,mirror,r,option,num_terms)
         [Z,coeffs] = FT_thermal_lens_from_coating_absorption(beam,mirror,r,num_terms);
         [c,p] = FT_approximate_paraboloid_for_thermal_lens(coeffs,mirror.X,mirror.a,beam.w);
         
-    % thermal lens from bulk absorption
+    % Thermal lens from bulk absorption
     elseif strcmp(option,'bulk')
        
         [Z,coeffs] = FT_thermal_lens_from_bulk_absorption(beam,mirror,r,num_terms);
diff --git a/src/new/thermal/FT_thermal_lens_from_bulk_absorption.m b/src/new/thermal/FT_thermal_lens_from_bulk_absorption.m
index d9527f571b8c9688689eea68cd3d006560890ae9..d157028148541d0782a7b10ab869d548f6c5c2e5 100644
--- a/src/new/thermal/FT_thermal_lens_from_bulk_absorption.m
+++ b/src/new/thermal/FT_thermal_lens_from_bulk_absorption.m
@@ -33,14 +33,13 @@ function [Z_bulk,coeffs] = FT_thermal_lens_from_bulk_absorption(beam,mirror,r,nu
     % Currently only valid for LG00
     p = 0;
     l = 0;  
-    
     a = mirror.a;
     h = mirror.h;
     K = mirror.K;
     dn_dT = mirror.dn_dT;
     Babs = mirror.Babs;
     X = mirror.X;
-    P = beam.P;
+    Power = beam.P;
     lambda = beam.lambda;
     w = beam.w;
     
@@ -48,21 +47,23 @@ function [Z_bulk,coeffs] = FT_thermal_lens_from_bulk_absorption(beam,mirror,r,nu
     r_abs = abs(r);
 
     % Constant
-    c = dn_dT * Babs*P*h / (pi*K); 
+    c = dn_dT * Babs*Power*h / (pi*K); 
     
     % 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);
+    zeta_n = FT_bessel_zeros(2*l,X,num_terms);
+    pn = FT_calculate_fourier_bessel_coefficients(p,l,X,w,a,num_terms);
     
     % Calculating sum of terms to approximate thermal lens
     Z_bulk=0;
     
     for s=1:num_terms
         
-        coeffs(s) = c * (1 - (2*X*a/(zeta_0s(s)*h) * sinh(zeta_0s(s)*h/(2*a))) / (zeta_0s(s)*sinh(zeta_0s(s)*h/(2*a)) + X*cosh(zeta_0s(s)*h/(2*a))));
-        coeffs(s) = coeffs(s) * exp(-zeta_0s(s)^2*w^2/(8*a^2)) / ((zeta_0s(s)^2+X^2) * (besselj(0,zeta_0s(s))^2));
+        gamma_n = zeta_n(s)*h/(2*a);
+        d1n = zeta_n(s)*sinh(gamma_n)+X*cosh(gamma_n);
         
-        Z_bulk = Z_bulk +  coeffs(s)*besselj(0,zeta_0s(s)*r_abs/a);
+        coeffs(s) = c * pn(s)/(zeta_n(s)^2) * (1-(X/gamma_n)*sin(gamma_n)/d1n);
+
+        Z_bulk = Z_bulk +  coeffs(s)*besselj(2*l,zeta_n(s)*r_abs/a);
     end
           
 
diff --git a/src/new/thermal/FT_thermal_lens_from_coating_absorption.m b/src/new/thermal/FT_thermal_lens_from_coating_absorption.m
index 93a4d5066f0716511de2b3a5914828dc2d70789f..734f2f8efd481afff0c7491ea475c7a41c2a0010 100644
--- a/src/new/thermal/FT_thermal_lens_from_coating_absorption.m
+++ b/src/new/thermal/FT_thermal_lens_from_coating_absorption.m
@@ -30,36 +30,41 @@
 function [Z_coat,coeffs] = FT_thermal_lens_from_coating_absorption(beam,mirror,r,num_terms)
     
     % Mirror and beam parameters 
-    % Currently only valid LG00
+    % Currently only valid LG00 
     p = 0;
-    l = 0; 
+    l = 0;
     a = mirror.a;
     h = mirror.h;
     K = mirror.K;
     dn_dT = mirror.dn_dT;
     Cabs = mirror.Cabs;
-    X = mirror.X
-    P = beam.P;
+    X = mirror.X;
+    Power = beam.P;
     lambda = beam.lambda;
     w = beam.w;
     
+    
     % Lens is axis-symmetric
     r_abs = abs(r);
 
     % Constant
-    c = dn_dT * Cabs*P / (pi*K); 
+    c = dn_dT * Cabs*Power / (pi*K); 
     
     % 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);
+    zeta_n = FT_bessel_zeros(2*l,X,num_terms);
+    pn = FT_calculate_fourier_bessel_coefficients(p,l,X,w,a,num_terms);
     
     % Calculating sum of terms to approximate thermal lens
     Z_coat=0;
     
     for s=1:num_terms
         
-        coeffs(s) = c*p_0s(s)*(1/zeta_0s(s)) * (1-exp(-zeta_0s(s)*h/a))/(zeta_0s(s)+X-(zeta_0s(s)-X)*exp(-zeta_0s(s)*h/a));
-        Z_coat = Z_coat +  coeffs(s)*besselj(0,zeta_0s(s)*r_abs/a);
+        % Coefficients
+        gamma_n = zeta_n(s)*h/(2*a);
+        d1n = zeta_n(s)*sinh(gamma_n)+X*cosh(gamma_n);
+        
+        coeffs(s) = c*pn(s)*(1/zeta_n(s)) * sinh(gamma_n)/d1n;
+        Z_coat = Z_coat +  coeffs(s)*besselj(2*l,zeta_n(s)*r_abs/a);
     end