addpath(genpath('/ligo/svncommon/SeiSVN/seismic/Common/MatlabTools'))
addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/


%% Actual ground motion

dt = 600;
start_time = 1127018405;
gps = start_time;
gwd = GWData();
gwd.make_kerberos_ready;

channels = {'H1:ISI-GND_STS_ETMY_Z_DQ,raw',...
    'H1:ISI-GND_STS_ETMX_Z_DQ,raw',...
    'H1:ISI-GND_STS_ITMY_Z_DQ,raw'};

gwd.site_info(2).server = 'nds2.ligo-wa.caltech.edu';
gwd.site_info(2).port = 31200;
 
[dataSTS, tSTS, info] = gwd.fetch(gps, dt, channels);

srate = 1/(tSTS(2) - tSTS(1));
nfft = srate * 8;




[DispEY.vel,DispEY.f] = pwelch(dataSTS(:,1),hanning(nfft),nfft/2,nfft,srate);
[DispEX.vel,DispEX.f] = pwelch(dataSTS(:,2),hanning(nfft),nfft/2,nfft,srate);
[DispIY.vel,DispIY.f] = pwelch(dataSTS(:,3),hanning(nfft),nfft/2,nfft,srate);



STScal = zpk([],0,1.6e-10);
STScalResp = abs(squeeze(freqresp(STScal,DispEY.f,'Hz')));


DispEY.disp = DispEY.vel .* STScalResp;
DispEX.disp = DispEX.vel .* STScalResp;
DispIY.disp = DispIY.vel .* STScalResp;


%%

G = 6.67e-11;    % m^3/(kg*s^2)
c = 250;         % m/s
rho_0 = 2500;    % kg/m^3
h = 1.5;         % m


a_NN_EY = 0.8 * 2*pi * G * rho_0 * DispEY.disp .* exp(-h * 2*pi * DispEY.f / c);
a_NN_EX = 0.8 * 2*pi * G * rho_0 * DispEX.disp .* exp(-h * 2*pi * DispEX.f / c);
a_NN_IY = 0.8 * 2*pi * G * rho_0 * DispIY.disp .* exp(-h * 2*pi * DispIY.f / c);

a_NN = 2 * a_NN_EY;  % Assuming all have same spectra, add incoherently
a_NN = sqrt(a_NN_EY.^2 + a_NN_EX.^2 + a_NN_IY.^2 + a_NN_IY.^2);  % IX and IY share an STS




Disp_NN = a_NN ./ (2*pi*DispEY.f).^2;

Strain_NN = Disp_NN / 4000;

%%

DARM = load('DARM.txt');  % from 23 Sept 2015




%%

loglog(DARM(:,1), DARM(:,2),...
    DispEY.f, Disp_NN,...
    'linewidth',2)
title('Estimated Newtonian noise, from measured Z-axis motion')
legend('DARM','Newtonian Noise')

xlabel('Frequency [Hz]')
ylabel('Displacement [m/\surdHz]')
set(gca,'FontSize',20)
axis([1, 500, 1e-22, 1e-15])

grid 
grid minor



