% base on LIGO-T0900209 equation 1

kB = 1.38e-23; % Boltzmann constant
T = 300; % temperature
Y = 7.27e10; % Young modulus, Quartz Glass for Optics Data and Properties, Heraeus data sheet
h = 35e-3; % [m] CP thickness VIR–0128A–12
n_s = 1.44963; % refraction index
r0 = 50e-3/sqrt(2); % [m] beam radius Braginsky definition
phi_s = 1/431; % CP loss angle, high due to rigid mounting, Effect of the CO2 scanning system on CP motion, Matteo Lorenzini and Paola Puppo
a = 350e-3/2; % CP radius


% frequency vector
F = logspace(1,3,1e3);
Omega = 2*pi*F;

% bessel function zeros
bessj1 = inline('besselj(1,x)');
for n = 1:10
  z(n) = fzero(bessj1,[n n+1]*pi);
end

% series expansion
k_m = z/a;

S_BR = 32*kB*T ./ (pi* Y*Omega) * ...
       sum( (cosh(k_m*h/2).^2) ./ (sinh(k_m*h) - k_m*h) .* exp(-k_m.^2*r0^2/2) ) * ...
       (n_s - 1)^2 * phi_s;

 sum( (cosh(k_m*h/2).^2) ./ (sinh(k_m*h) - k_m*h) .* exp(-k_m.^2*r0^2/2) )

% plot result
figure;
loglog(F, sqrt(S_BR));
grid;
xlabel('Frequency (Hz)')
ylabel('Optical length (m/rtHz)')

% include arm length, effect of arm cavity finesse, and the fact that there are two CP
figure;
loglog(F, sqrt(S_BR)/3e3*5e-3*sqrt(2));
grid;
xlabel('Frequency (Hz)')
ylabel('Strain (1/rtHz)')