% Vertical Blade Spring Design
% which has a dual-triagular shape
% see Stefan Gossler's Thesis, as
% well as LIGO T030107, and Calum and 
% Norna's blade spreadsheet
%
% Bram Slagmolen - July 09
%

fname = ['TT-' date '.txt'];
delete(fname);
diary(fname);   % turn on writing command window output saving to file

disp(' ');

material = 'BeCu';
%material = 'SS';

t = 0.254e-3;     % thickness of the Blade Spring, 306e-6 m, 0.006 inch  (0.152 mm), 0.01 inch (0.254 mm) or 0.016 inch (0.406 mm)
a = 26e-3;    % width of the base of the Blade Spring
a1 = 1e-10; % width of the tip of the Blade, the 1e-10 constant is to keep Beta NOT 1.
L = 50e-3; % length of the Leaf Spring

Beta = a1/a;
alpha = ((3/(2*(1-Beta))) * (3 - (2/(1-Beta))*(1+Beta^2 * log10(Beta)/(1-Beta)))); % GEO value 1.36 - 1.45
alpha = 1.5;      % this is the shape factor, which for a triangular blade is 1.5 and 
                % a rectangle 1.
g = 9.81;

%if material == 'ALU'
    E_al = 69e9;    % Youngs modulus of aluminium
    Y_al = 400e6;   % Yield strength
%elseif material == 'SS'
    E_ss = 210e9;   % Youngs modulus of Steel
    Y_ss = 600e6;   % Yield strength of Stainless Steel (approx.)
%elseif material == 'BeCu'    
    E_becu = 131e9; % Youngs molulus Berrylium-Copper
    Y_becu = 620e6; % Yield Strength Berrylium-Copper
%else material == 'Maraging'
    E_maraging = 186e9; % Youngs modulus for Marval18 Maraging Steel
    Y_maraging = 186e9; % Yield Strength o Maraging Steel
%end

Y_pw = 2200e6;  % Yield strength of piano wire
Rho_becu = 8100; % density of berillium copper

II = (a*t^3)/12; % second moment of inertia of the leaf spring

disp(fname);
disp('----------------------------');
disp(['blade material: ', material]);
disp(['blade thickness: ', num2str(t*1e3), ' mm (',num2str(t*1e3/25.4),' inch)']);
disp(['width at the base: ', num2str(a*1e3), ' mm (',num2str(a*1e3/25.4),' inch)']);
disp(['width at the tip: ', num2str(a1*1e3), ' mm (',num2str(a1*1e3/25.4),' inch)']);
disp(['blade length: ', num2str(L*1e3), ' mm (',num2str(L*1e3/25.4),' inch)']);
disp(['Beta (width at tip/width at base): ', num2str(Beta)]);
disp(['Alpha (1.3 - 1.45): ', num2str(alpha)]);
disp(['second moment of inertia (cross section): ', num2str(II), ' m^4']);
disp(' ');

% for a simple blade system
E = E_becu;
Y = Y_becu;

% Effective spring constant
k = (E* a * t^3)/(4*L^3 * alpha); % stiffness [N/m]
disp(['effective sping constant: ', num2str(k), ' N/m']);

% Effective suspended mass
mirrormass = 0.124;
clampmass = 0.6e-3;
m = mirrormass/2 + clampmass;
disp(['effective mass: ', num2str(m), ' kg']);

% Effective resonance frequency
f0 = sqrt(E*a*t^3 / (16*pi^2 * m * L^3 * alpha)); % from T030107
disp(['effective resonance frequency: ', num2str(f0), ' Hz']);

% Blade tip deflection
d = 4*m*g*L^3 * alpha / (E*a*t^3);  % from T030107
disp(['deflection: ', num2str(d*1e3), ' mm']);

% maximum stress
Tm = 6 * g * m * L / (a * t^2);
disp(['maximum stress: ', num2str(Tm/1e6), ' MPa']);


% Rc = E * a * t^3 / (12 * m * g * L); % from T030107
% From Norna's spreadsheet -> d = Rc * [1 - cos(L/R)], 
% -> shoudln't this be d = Rc * [1 - cos(L/2R)]?
% using taylor series cos(x)=1 - x^2/2! + x^4/4!
% Rc = | L^2 / (2*d) |
Rc = L^2 / 2 / d;
disp(['pre-curved radius at rest: ', num2str(Rc), ' m']);

% blade base to tip length under load
% l = Rc * sqrt( 2 - 2*cos(L/Rc) )
l = Rc * sqrt(2 - 2*cos(L/Rc)  );
disp(['distance base to tip underload: ',num2str(l*1e3),' mm']);

% angle of the blade tip in respect to the base
theta = 90 - acos(L/Rc)*180/pi;
disp(['angle at blade tip is: ', num2str(theta), ' deg']);

% Angle of the suspension wire from the suspension point (blade tip)
% to the optic barrel. The d-yaw at suspension point is 35mm
% while the optic diamter is 3", and the supsension length is 140mm
dyaw = 35e-3;
doptic = 3*25.4e-3;
Lpend = 140e-3;
theta2 = ((doptic-dyaw)/2 /Lpend)*180/pi; 
disp(['suspension wire angle: ', num2str(theta2), ' deg']);

% Angle which the blade will be mounted so the suspension wire
% leaves the blade at 90 deg (or as close as possible).
theta3 = theta - theta2;
disp(['angle blade base is mounted (90 deg of wire to clamp): ', num2str(theta3), ' deg']);


diary off;  % turn off diary, e.g. saving output to file

% Using Dennis calculations
B1 = a/a1;

lk0 = 1.3;
f00 = (1/2/pi) * (t/L^2) * lk0 * sqrt(E*g/12/Rho_becu);


