% Calculation of suspension thermal noise.
%
% f:                    frequency vector
% par:                  Speedmeter parameter set
% mirror:               name of mirror to calculate for, e.g. 'M1a'
% beamOffset:           offset of laser beam from centre of mirror
% vertToHorCoupling:    coupling from vertical motion to horizontal motion
% nViolin:              number of violin modes to calculate
%
% Script originally by Andrew Spencer, 2015-03-05. The calculations for
% suspension thermal noise are all from his original script, on the labbok
% at https://arran.physics.gla.ac.uk/wp/speedmeter/2015/04/14/final-project-report/
%
% Made compatible with the Speedmeter parameter set by Sean Leavey,
% 2015-11-20.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Change Log %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%% 09 Feb 2016, JSH
    % Added missing factor rho to formula for tau
    % Added square in formula for TE loss
    % Added missing dot in formula for TE loss
    % Substituted ss with 62.45e6 (sigma in that formula is the static
    %   stress in one fibre and not the poisson ratio)
    % Added /2 in formula for area moment of inertia
    
    
function totalNoise = calcSusThermNoise(f, par, mirror, beamOffset, ...
    vertToHorCoupling, nViolin)
    %% Paths
    
    % this script's directory
    directory = fileparts(mfilename('fullpath'));
    
    % validation functions
    addpath(genpath([directory, '/../../IfoSimulations/Libs/Validation/']));

    %% Input checking
    
    if ~isFrequencyVector(f)
        error('Specified frequency vector is invalid');
    end
    
    if ~isParamSet(par)
        error('Specified parameter set is invalid');
    end
    
    if ~ischar(mirror)
        error('Specified mirror must be a string');
    end
    
    %% Physical constants

    % gravity
    g = par.constants.g;
    
    % temperature
    T = par.thermal.ETcryo;
    
    % Boltzmann constant
    kB = par.constants.kB;

    %% Mirror properties
    
    % get mirror properties from parameter set
    if par.Mirror.isKey(mirror)
        optic = par.Mirror(mirror);
    elseif par.BS.isKey(mirror)
        optic = par.BS(mirror);
    else
        error('Specified mirror is not part of the specified parameter set');
    end
    
    % check that fibre geometry exists for this optic
    if ~isfield(optic, 'fibregeometry')
        error(['Specified mirror does not contain fibre geometry information.', ...
            'Add this to the parameter set if you wish to calculate noise for this mirror.']);
    end
    
    % fibre properties
    numFibres = optic.fibregeometry.nFibres;
    fibLength = optic.fibregeometry.fibreLength;
    fibDiameter = optic.fibregeometry.fibreDiameter;
    %fibDiameter = 10e-6;
    
    % fibre attachment points
    attVertDist = optic.fibregeometry.attVertDist;
    attHorDist = optic.fibregeometry.attHorDist;
    attTransDist = optic.fibregeometry.attTransDist;
    
    % mirror geometry
    mirrorMass = optic.mass
    mirrorRadius = optic.diameter / 2;
    mirrorThickness = optic.thickness;
    
    % material constants
    phi = optic.bulkmaterial.phi;
    Y = optic.bulkmaterial.Y;
    rho = optic.bulkmaterial.rho;
    C = optic.bulkmaterial.C;
    alpha = optic.bulkmaterial.alpha;
    Beta = optic.bulkmaterial.beta;
    ss = optic.bulkmaterial.sigma;
    K = optic.bulkmaterial.k;
    h = optic.bulkmaterial.dissdepth;

    %% Moment on inertia through central axis
    Iy = 0.25*mirrorMass*(mirrorRadius^2)+ 0.0833*mirrorMass*(mirrorThickness^2);
    Ix = 0.5*mirrorMass*(mirrorRadius^2);

    %% Flexure length
    a = sqrt((Y*Ix)/(mirrorMass*g));

    %% Spring constant of one fiber
    k = (Y*pi*(fibDiameter/2)^2)/(fibLength);
    Y
    fibLength
    fibDiameter
    %% Matrix elements
    k11 = -((mirrorMass*g*attVertDist)/(Iy))-((mirrorMass*g*(attVertDist^2))/(Iy*fibLength))-((numFibres*k*(attHorDist^2))/(Iy));
    k12 = (mirrorMass*g*attVertDist)/(Iy*fibLength);
    k21 = (g*attVertDist)/fibLength;
    k22 = -g/fibLength;
    A = [k11, k12;k21, k22];

    s11 = -((mirrorMass*g*attVertDist)/(Ix))-((mirrorMass*g*(attVertDist^2))/(Ix*fibLength))-((numFibres*k*(attTransDist^2))/(Ix));
    s12 = (mirrorMass*g*attVertDist)/(Ix*fibLength);
    s21 = -(g*attVertDist)/fibLength;
    s22 = -g/fibLength;
    B = [s11, s12;s21, s22];

    %% Resonant frequencies of the different modes of motion
    flongpitch = (sqrt(abs(eig(A))))./(2*pi);
    fsideroll = (sqrt(abs(eig(B))))./(2*pi);
    fvert = (sqrt((numFibres*k)/(mirrorMass)))/(2*pi);
    fyaw = (sqrt((mirrorMass*g*(attTransDist^2+attHorDist^2))/(fibLength)))/(2*pi);

    % assign the eigenfrequencies to the correct mode
    ap = sqrt(attVertDist^2 + attHorDist^2);

    if ap == 0
        flong =flongpitch(1);
        fpitch = flongpitch(2);
    else
        fpitch = flongpitch(1);
        flong = flongpitch(2);
    end

    ap2 = sqrt(attVertDist^2 + attTransDist^2);

    if ap2 == 0
        fside = fsideroll(1);
        froll = fsideroll(2);
    else
        froll = fsideroll(1);
        fside = fsideroll(2);
    end

    % point mass pendulum mode
        f0 = (1/(2*pi))*(sqrt(g/fibLength));

    %% Precalculations

    % cross-sectional moment of inertia
    I = (pi*(fibDiameter/2)^4)/4

    % tension in each fiber
    t = (mirrorMass*g)/numFibres

    % dissipation dilution
    D = (1/numFibres/2)*((2*mirrorMass*g*fibLength)/(sqrt(t*Y*I)))

    % weld loss;
    %phi_w = 5.8e-7;     % Value as an estimate from aLIGO fibre
    phi_w = 0;
    % mech loss of epoxy
    phi_ep = 0.1; 

    % thermoelastic loss
    tau = (C*(fibDiameter^2)*rho)/(4.32*pi*K);
    F = 2.*pi.*f.*tau;
    phi_te = ((Y.*T)./(rho.*C)).*(alpha-(62.45e6*Beta)./Y).^2.*((F)./(1+(F.^2)));

    % surface loss
    %phi_s = (8*h*phi)/fibDiameter;
    phi_s = (8*5e-13)/fibDiameter;
    % total mechanical loss
    phi_tot = (phi+phi_w+phi_te+phi_s)./D;

    % prefactor
    c = (kB*T)/(2*(pi^3));

    % isplacement due to vertical mode
    TNvert2 = (c./(mirrorMass.*f)).*((phi_tot.*fvert.^2)./((phi_tot.^2).*(fvert.^4)+((fvert.^2)-(f.^2)).^2));

    % vertical mode contribution to longitudinal laser noise in
    % amplitude spectral density
    TNvert = vertToHorCoupling.*(sqrt(TNvert2));

    % displacement due to londitudinal and pitch modes
    TNangpitch2 = (c./(mirrorMass.*f)).*((phi_tot.*fpitch.^2)./((phi_tot.^2).*(fpitch.^4)+((fpitch.^2)-(f.^2)).^2));
    TNangpitch = sqrt(TNangpitch2);
    TNlong2 = (c./(mirrorMass.*f)).*((phi_tot.*flong.^2)./((phi_tot.^2).*(flong.^4)+((flong.^2)-(f.^2)).^2));

    % convert angular displacment noise to linear amplitude spectral density displacement noise for
    % pitch mode
    TNpitch = (beamOffset.*(TNangpitch))+((mirrorThickness.*(TNangpitch2))./4);
    TNlong = sqrt(TNlong2);

    % yaw mode
    TNyaw2 = (c./(mirrorMass.*f)).*((phi_tot.*fyaw.^2)./((phi_tot.^2).*(fyaw.^4)+((fyaw.^2)-(f.^2)).^2));
    TNyaw = (beamOffset.*(sqrt(TNyaw2)))+((mirrorThickness.*(TNyaw2))./4);

    % sideways mode
    TNside2 = (c./(mirrorMass.*f)).*((phi_tot.*fside.^2)./((phi_tot.^2).*(fside.^4)+((fside.^2)-(f.^2)).^2));
    TNside = vertToHorCoupling.*(sqrt(TNside2));
    TNroll2 = (c./(mirrorMass.*f)).*((phi_tot.*froll.^2)./((phi_tot.^2).*(froll.^4)+((froll.^2)-(f.^2)).^2));
    TNroll = vertToHorCoupling.*(sqrt(TNroll2));

    % empty matrix whose nth row will be filled with the lineshape of
    % nth violin mode
    TN = zeros(nViolin, length(f));

    for n = (1:1:nViolin);
        % define resonant frequency and reduced mass of the nth violin mode
        f0v = (n./(fibLength.*fibDiameter)).*(sqrt(t./(pi.*rho)));
        mu = 0.5*mirrorMass*((f0v/f0)^2);
        TNv2 = 4.*(c./(mu.*f)).*((phi_tot.*f0v.^2)./((phi_tot.^2).*(f0v.^4)+((f0v.^2)-(f.^2)).^2));
        TNv = sqrt(TNv2);
        TN(n,:) = TNv;
    end

    TNviolin = sum(TN);

    % total noise
    totalNoise = TNviolin + TNlong + TNpitch + TNvert + TNyaw + TNside;
end