%%  function to calculate total suspension thermal noise depending on frequency
%   calling frequency, constants from list_constants.m


%% funciton total_sus_therm

function [total_sus_therm] = total_sus_therm(freq)

list_constants;

T = T_mir;            % temperature
m = M4;               % mass suspension
L = F_l;              % length suspension
rad1 = L4_rad1;       %for 40kg      %radius of thick section of fibre (to cancel thermoelastic)
rad2 = L4_rad2;       %for 40kg      %radius of thin section of fibre

Y = constants.Y;
rho = constants.rho;
C = constants.C;
kB = constants.kB;

kappa = constants.kappa;             % thermal conductivity
alpha = constants.alpha;             % thermal expansion coefficient silica
beta = constants.beta;               % Young's modulus variation of silica (1/Y dY/dT)

Iarea = pi.*rad1^4/4;                % second moment of areafor bending region
hphi_s=6.15e-12;                     % surface loss
phi_weld=5.8e-7;                     % weld loss
weld_energy=7e-2;                    % weld energy
tension=(m/4).*9.81;                 % tension
kl=m.*9.81/L;                        % longitudinal spring constant
kv=Y.*pi.*(rad2)^2/L;                % vertical spring constant
wl=sqrt(9.81/L);                     % longitudinal w
wv=sqrt(2).*sqrt(kv/(m/4));          % vertical w (sqrt(2) higher as both masses free, assumes masses identical=> need to change 


% violin mode (just fundamental mode)

% mass per unit length
mu=pi.*rad2^2.*rho;
mass_fibre=mu.*L;

speed=sqrt(tension/mu);
fviolin1=(speed/(2.*L));
wviolin1=2.*pi.*fviolin1;

%factor=0.5;     % half dilution for bending at both ends
%FEA=1/3;      % factor to account for variation between analytical and FEA (about 1.8)

%dilution_l=2.*L.*sqrt(tension/(Y.*Iarea)).*factor.*FEA     % ideal dilution for longitudinal
%dilution for L4=1.2m
if val==0
    dilution_l=91
else
if m==40
    dilution_l=269.1269*(L/1.2) %40kg
    %dilution_l=91;
elseif m==60
    dilution_l=217.7803*(L/1.2) %60kg
elseif m==80
    dilution_l=187.5951*(L/1.2) %80kg
elseif m==100
    dilution_l=166.6399*(L/1.2) %100kg
elseif m==120
    dilution_l=151.736*(L/1.2)  %120kg
elseif m==140
    dilution_l=140.3183*(L/1.2) %140kg
elseif m==160
    dilution_l=139.6907*(L/1.2) %160kg
elseif m==180
    dilution_l=123.7676*(L/1.2) %180kg
elseif m==200
    dilution_l=117.6563*(L/1.2) %200kg
end
end

cross_coupling=1/1000;                                      % "dilution" for vertical => cross coupling of 0.1%
dilution_violin=0.5.*L.*sqrt(tension/(Y.*Iarea))           % ideal dilution for violin (approximate solution, needs checking)

% define frequency points
freq=frequency;     % freq = logspace(log10(fmin),log10(fmax),N); would be another interesting option
omega=2.*pi.*freq;

% assume idealistic loss and that all energy is stored in the 800um region
phi_surf=(8.*hphi_s/(2.*rad1));         % phi surface
phi_weld=(phi_weld).*weld_energy;       % phi weld

% phi thermoelastic
tau1=1/(4.32.*pi).*rho.*C.*rad1^2/kappa;
sigma1=tension/(pi.*(rad1)^2); 
val1=(Y.*T)/(rho.*C).*(alpha-sigma1.*(beta/Y))^2;
%--------------------------------------------------------------------

n=length(freq)
for j=1:n
% phi thermoelastic
phi_thermo=(omega(j)*tau1)/(1+(omega(j)*tau1)^2)*val1;
% phi bulk
phi_bulk=4.1e-12*freq(j)^0.77;

% longitudinal loss
phi_term_l=(phi_bulk+phi_thermo+phi_surf+phi_weld)/dilution_l;
% vertical loss (no thermoelastic, surface/2)
phi_term_v=(phi_bulk+phi_surf/2);
% violin loss (2* undiluted longitudinal)
phi_term_violin1=2*(phi_term_l*dilution_l)/dilution_violin;
    
phi_total_l(j)=phi_term_l;   
phi_total_v(j)=phi_term_v;   
phi_total_violin1(j)=phi_term_violin1;   
phi_total_surf(j)=phi_surf;
phi_total_bulk(j)=phi_bulk;
phi_total_thermo(j)=phi_thermo;
phi_total_weld(j)=phi_weld;

% calculate longitudinal thermal displacement noise
xthermal(j)=sqrt((4*kB*T*phi_total_l(j)*wl^2)/((m*omega(j))*((wl^2-omega(j)^2)^2+(phi_total_l(j)*wl^2)^2)));
% calculate vertical thermal displacement noise
zthermal(j)=cross_coupling*sqrt((4*kB*T*phi_total_v(j)*wv^2)/((m*omega(j))*((wv^2-omega(j)^2)^2+(phi_total_v(j)*wv^2)^2)));
% calculate fundamental violin mode thermal displacement noise
v1_thermal(j)=sqrt((8*kB*T*mass_fibre*phi_total_violin1(j)*wviolin1^2)/((pi^2*m^2*omega(j))*((wviolin1^2-omega(j)^2)^2+(phi_total_violin1(j)*wviolin1^2)^2)));
total(j)=sqrt(xthermal(j)^2+zthermal(j)^2+v1_thermal(j)^2);
total_sus_therm(j) = total(j);
x_therm(j) = xthermal(j);
z_therm(j) = zthermal(j);
v_therm(j) = v1_thermal(j);
end