%% Function to calculate longitudinal seismic noise (approx)
%% Function long_seis_approx = long_seis_approx(freq)
function [long_seis_approx] = long_seis_approx(freq)

list_constants;

%% Initializations

total_wire_length = total_length;      % length of suspension in m
L4 = F_l;                      % fiber length in m; effective value, should change these to CM to CM value
L1 = L_1; % top wire length; effective value, should change these to CM to CM value
L2 = L_2;                       % top mass to UIM wire length; effective value, should change these to CM to CM value
L3 = L_3;                       % UIM to PUM wire length; effective value, should change these to CM to CM value


m1 = M1; % top mass mass in kg
m2 = M2; % top mass mass in kg
m3 = M3; % top mass mass in kg
m4 = M4;   % top mass mass in kg
g  = constants.g;  % gravitational acceleration in m/s^2
Y3 = constants.Y; % fiber modulus of elastcity; modified by Madeleine Waller with info from Eric Gustafson 8/13/12 (also consistent with 'Temperature dependence of the force sensitivity of silicon cantilevers', which gives 166 GPa at 295 K)                       
fiber_stress = 0.768881e+09; % assumes a fiber area of 0.25 mm^2 with a 143 kg test mass. This area is halfway between Alastair's 0.4 mm^2 and Rana's 0.1 mm^2
freq=frequency;
%freqval=10; % Hz, frequency to evaluate sus longitudinal isolation
aLIGO_isolation_factor = 1.0875e-7; % m/m between the suspension point and test mass along longitudinal axis

% Longitudinal Isolation with initial parameters

n=length(freq);
for j=1:n
xg_approx(j) = ((g.^4)./((2*pi*freq(j)).^8))*(1./(L1*L2*L3*L4))*(m1+m2+m3+m4)*(m2+m3+m4)*(m3+m4)./(m1*m2*m3); % approximate longitdudinal ilsolation from TF asymptote, equation 1 in T1300786-v7.
long_seis_approx(j) = xg_approx(j);

end