printFigs = true;
fetchdata = true;
sampf = 8; % sampling frequency in Hz

% Define Channels
if fetchdata
% conn = nds2.connection('nds2.ligo-wa.caltech.edu', 31200);
Channels = {'H1:ISI-GND_STS_ETMX_X_DQ',...
'H1:ISI-GND_STS_ETMX_Y_DQ',...
'H1:ISI-GND_STS_ETMX_Z_DQ',...
'H1:ISI-GND_STS_ETMY_X_DQ',...
'H1:ISI-GND_STS_ETMY_Y_DQ',...
'H1:ISI-GND_STS_ETMY_Z_DQ',...
'H1:ISI-GND_STS_ITMY_X_DQ',...
'H1:ISI-GND_STS_ITMY_Y_DQ',...
'H1:ISI-GND_STS_ITMY_Z_DQ',...
'H1:ISI-GND_BRS_ETMY_RX_IN1_DQ',...
};

%Starttime = 1097105000; % windy data
%Starttime = 1097560Rot_000;
%Starttime=1143572417; %April 1, 2016 12:00:00
Starttime=1143962787;
%Starttime = 1143824000; % Windy Data ~25 mph 04/04/2016
duration = 3000;

% numfetches=length(Channels);
% rawdata=zeros(duration*sampf,numfetches);
% %Fetch data
% for i=1:numfetches
% buffers = conn.fetch(Starttime, Starttime+duration, Channels(i));
% %Down sample to desired value (presently 8 Hz)
% rawdata(:,i) = decimate(double(buffers(1).getData), buffers(1).getChannel.getSampleRate/sampf);
% conn.clearCache;
% end
rawdata = get_data(Channels,'raw',Starttime,duration);
rawETMXX = decimate(rawdata(1,1).data, rawdata(1,1).rate/sampf);
rawETMXY = decimate(rawdata(1,2).data, rawdata(1,2).rate/sampf);
rawETMXZ = decimate(rawdata(1,3).data, rawdata(1,3).rate/sampf);
rawETMYX = decimate(rawdata(1,4).data, rawdata(1,4).rate/sampf);
rawETMYY = decimate(rawdata(1,5).data, rawdata(1,5).rate/sampf);
rawETMYZ = decimate(rawdata(1,6).data, rawdata(1,6).rate/sampf);
rawITMYX = decimate(rawdata(1,7).data, rawdata(1,7).rate/sampf);
rawITMYY = decimate(rawdata(1,8).data, rawdata(1,8).rate/sampf);
rawITMYZ = decimate(rawdata(1,9).data, rawdata(1,9).rate/sampf);
rawBRS = decimate(rawdata(1,10).data, rawdata(1,10).rate/sampf);
mydata=[rawETMXX rawETMXY rawETMXZ rawETMYX rawETMYY rawETMYZ rawITMYX rawITMYY rawITMYZ rawBRS];
save('GPS1143962787_6_9Earthquake.mat','mydata')
end
if exist('GPS1143962787_6_9Earthquake.mat','file')
myfile = load('GPS1143962787_6_9Earthquake.mat');
mydata = myfile.mydata;
rawETMXX= mydata(:,1);
rawETMXY = mydata(:,2);
rawETMYZ = mydata(:,3);
rawETMYX = mydata(:,4);
rawETMYY = mydata(:,5);
rawETMYZ = mydata(:,6);
rawITMYX = mydata(:,7);
rawITMYY = mydata(:,8);
rawITMYZ = mydata(:,9);
rawBRS= mydata(:,10);
end

set(0,'DefaultAxesFontSize',20);

Sttime = 001 * sampf;
Totaltime = 2990; %duraiton in sec
Endtime = Sttime + Totaltime * sampf;
localg = 9.8;
BRSscale=.85;

ETMXX=1e-9 *rawETMXX(Sttime:Endtime);
ETMXY=1e-9 *rawETMXY(Sttime:Endtime);
ETMXZ=1e-9 *rawETMXZ(Sttime:Endtime);
ETMYX=1e-9 *rawETMYX(Sttime:Endtime);
ETMYY=1e-9 *rawETMYY(Sttime:Endtime);
ETMYZ=1e-9 *rawETMYZ(Sttime:Endtime);
ITMYX=1e-9 *rawITMYX(Sttime:Endtime);
ITMYY=1e-9 *rawITMYY(Sttime:Endtime);
ITMYZ=1e-9 *rawITMYZ(Sttime:Endtime);
BRS=1e-9* BRSscale*rawBRS(Sttime:Endtime);

Navg = 5;

%% BRS parameters
Mtotal = 4.5; % Total mass in kg
Ibar = 0.61; % Moment of Inertia in kg m^2
localg = 9.8; % local gravitation acceleration
doffset = .5e-6; % CoM offset from pivot in m. d is plus if CoM below pivot
resfreq = 2*pi*7.7e-3;
Qbar = 3000;
Rot_time = transpose(1/sampf * (0:1:length(BRS)-1));

%% torque computation
[bb,aa] = butter(4,[2*0.03/sampf 2*0.100/sampf],'bandpass');
% %T240 response inversion filter
T240InvertFilt = zpk(-2*pi*[pairQ(8.2e-3,0.7)],-2*pi*[0 0],1);
T240InvertFilt = 1*T240InvertFilt/abs(freqresp(T240InvertFilt,2*pi*100));
%
% %BRS response inversion filter
BRSInvertFilt = zpk(-2*pi*[pairQ(7.7e-3,4000)],-2*pi*[0 0],1);
BRSInvertFilt = 1*BRSInvertFilt/abs(freqresp(BRSInvertFilt,2*pi*100));
%
% %Filters to differentiate and integrate
DiffFilt = zpk(-2*pi*[0], -2*pi*2,1);
DiffFilt = 1*DiffFilt/abs(freqresp(DiffFilt,2*pi*0.1592));
%
IntFilt = zpk(-2*pi*[], -2*pi*5e-4,1);
IntFilt = 1*IntFilt/abs(freqresp(IntFilt,2*pi*0.1592));
%
% %Apply filters
T240cal_vel = lsim(T240InvertFilt,ETMYZ,Rot_time);
%
T240cal_acc = lsim(DiffFilt,ETMYY,Rot_time);
%
T240cal_disp = lsim(IntFilt,T240cal_vel,Rot_time);
%
BRScal_out = lsim(BRSInvertFilt,BRS, Rot_time);
BRS_out=filter(bb,aa,BRScal_out);

ETMYZ_out=filter(bb,aa,T240cal_disp);
%
% ResidualAcc = T240cal_acc + BRScal_out + Mtotal*localg/Ibar*doffset*T240cal_disp;

%% Calculate ASD

[AETMXX, ~] = asd2(ETMXX,1/sampf, Navg, 1, @hann);
[AETMXY, ~] = asd2(ETMXY,1/sampf, Navg, 1, @hann);
[AETMXZ, ~] = asd2(ETMXZ,1/sampf, Navg, 1, @hann);
[AETMYX, ~] = asd2(ETMYX,1/sampf, Navg, 1, @hann);
[AETMYY, ~] = asd2(ETMYY,1/sampf, Navg, 1, @hann);
[AETMYZ, ~] = asd2(ETMYZ,1/sampf, Navg, 1, @hann);
[AITMYX, ~] = asd2(ITMYX,1/sampf, Navg, 1, @hann);
[AITMYY, ~] = asd2(ITMYY,1/sampf, Navg, 1, @hann);
[AT240Yaccel_cal, ~] = asd2(T240cal_acc,1/sampf, Navg, 1, @hann);
[AITMYZ, ~] = asd2(ITMYZ,1/sampf, Navg, 1, @hann);
[ABRS, ~] = asd2(BRS,1/sampf, Navg, 1, @hann);
[AETMYZ_out, Af] = asd2(ETMYZ_out,1/sampf, Navg, 1, @hann);

% [ADiffY, ~] = asd2((ETMY-ITMY),1/sampf, Navg, 1, @hann);
% [ACommY, ~] = asd2((ETMY+ITMY)/2,1/sampf, Navg, 1, @hann);
% [ADiffX, ~] = asd2((ETMX-ITMX),1/sampf, Navg, 1, @hann);
% [ACommX, frot] = asd2((ETMX+ITMX)/2,1/sampf, Navg, 1, @hann);
% [RY,ADiffOY,~,Tyy] = sccs2(ETMY,ITMY,1/sampf, Navg, 1, @hann);
% [RX,ADiffOX,fO,Txx] = sccs2(ETMX,ITMX,1/sampf, Navg, 1, @hann);

%% Calculate Coherence
Navg2 = 9;

[CohETMYX_BRS, ~] = coh2(ETMYX,BRS,1/sampf, Navg2, 1, @hann);
[CohETMYY_BRS, ~] = coh2(ETMYY,BRS,1/sampf, Navg2, 1, @hann);
[CohETMYZ_BRS, Cf] = coh2(ETMYZ,BRS,1/sampf, Navg2, 1, @hann);

%% Transfer Function
Navg3 = 17;
%[Tyy,~] = tfe2(ETMYZ_out, 1.3520e+04*BRS_out,1/sampf, Navg3);
[Tyy,~] = tfe2(ETMYZ,ITMYZ ,1/sampf, Navg3);
[Txx,F] = tfe2(ETMXZ,ITMYZ,1/sampf, Navg3);
figure(5)
l=semilogx(F,abs(Tyy),F,abs(Txx));
set(l,'LineWidth',2);
set(gca,'FontSize',18)
xlim([3e-3 .2])
xlabel('frequency (Hz)');
ylabel('Magnitude');
legend('ETMY_Z and ITMY_Z','ETMX_Z and ITMY_Z');
grid on

%% Phase
figure(6)
l=semilogx(F,angle(Tyy)*180/pi,F,angle(Txx)*180/pi);
set(l,'LineWidth',2);
set(gca,'FontSize',18)
xlim([3e-3 4])
xlabel('frequency (Hz)');
ylabel('Phase');
legend('ETMY_Z and BRS','ETMX_Z and ETMY_Z');
grid on
%% Plots

figure(1)
ll=loglog(Af,AETMXX,Af,AETMXY,Af,AETMXZ);
% ll=loglog(frot, ABRSdrift, frot, ABRSout,frot, AT240YV.*(2*pi*frot'),...
% frot,AT240Yaccel_cal.*Mtotal*localg/Ibar*abs(doffset)./(2*pi*frot').^2);
set(ll,'LineWidth',2);
set(gca,'FontSize',18)
xlim([3e-3 4])
% ylim([1e-10 1e-6])
%ylim([1e-10 1e-6])
set(gca,'YTick',10.^(-10:-4))
xlabel('frequency (Hz)');
ylabel('ASD (m/s/rt(Hz))');
% ylabel('ASD (rad/rt(Hz))');
title('STS');
% legend('BRS-2 rX in (raw)','BRS-2 rX OUT','STS Y (raw) * w^2/g','Acceleration coupling with d = 1e-6 m',...
% 'Location','northeast'); %'IPS X,RY noise')
% legend('DiffY','DiffX','DiffY Optimal','DiffX Optimal','CommX','CommY');
legend('ETMX_X','ETMX_Y','ETMX_Z')
grid on

figure(2)
ll=loglog(Af,AETMYZ,Af,ABRS,Af,AT240Yaccel_cal.*Mtotal*localg/Ibar*abs(doffset)./(2*pi*Af').^2);
set(ll,'LineWidth',2);
set(gca,'FontSize',18)
xlim([3e-3 4])
ylim([1e-11 1e-2])
set(gca,'YTick',10.^(-10:-2))
xlabel('frequency (Hz)');
ylabel('ASD (m/s/rt(Hz) or rad/rt(Hz))');
title('STS');
legend('STS_Z ETMY','BRS ETMY','BRS possible acceleration coupling')
grid on

figure(3)
ll=loglog(Af,AITMYX,Af,AITMYY,Af,AITMYZ);
set(ll,'LineWidth',2);
set(gca,'FontSize',18)
xlim([3e-3 4])
set(gca,'YTick',10.^(-10:-4))
xlabel('frequency (Hz)');
ylabel('ASD (m/s/rt(Hz))');
title('STS');
legend('ITMY_X','ITMY_Y','ITMY_Z')
grid on

figure(8)
ll=loglog(Af,ABRS* 1.3520e+04 ,Af,AETMYZ_out);
set(ll,'LineWidth',2);
set(gca,'FontSize',18)
ylim([1e-9 1e-3])
xlim([5e-3 5e-1])
set(gca,'YTick',10.^(-11:-3))
xlabel('frequency (Hz)');
ylabel('ASD (m/rt(Hz))');
title('BRS');
legend('BRS','ETMYZ')
grid on

figure(4)
%subplot(212)
lll = semilogx(Cf, CohETMYX_BRS,Cf,CohETMYY_BRS,Cf,CohETMYZ_BRS); %, frot2, CohBRSoutT240);
set(lll,'LineWidth',2);
xlim([3e-4 4])
xlabel('frequency (Hz)');
ylabel('coherence');
%legend('BRS and PEM tilt','T240 and PEM tilt','BRSout and T240X','Location','southeast');
legend('ETMY_X and BRS','ETMY_Y and BRS','ETMY_Z and BRS');
grid on

tim=(1:length(ETMYZ_out))/8;

figure(7)
plot(tim,ETMYZ_out+5e-5,tim,BRS_out*13250);
hold off
xlim([1.2e2 3e3])
ylim([-3e-5 1e-4])
legend('STS_Z ETMY','BRS ETMY')
ylabel('Displacement (m)')
xlabel('Time (s)')
grid on
%% seismometer plot
% figure(3)
% lllll = loglog(frot, AT240XV./(2*pi*frot'),frot, AT240YV./(2*pi*frot'),...
% frot, AT240ZV./(2*pi*frot')); %frot,Ibar*ABRSsubaccel);
% %frot,Ibar*ABRSsubaccel);
% set(lllll,'LineWidth',2);
% set(gca,'FontSize',18)
% xlim([3e-4 1])
% ylim([1e-10 1e-3])
% set(gca,'YTick',10.^(-10:-3))
% xlabel('frequency (Hz)');
% ylabel('velocity ASD (m/s/rt(Hz))');
% title('T240 seismometer');
% legend('X','Y','Z',...
% 'Location','northeast')%'IPS X,RY noise')
% grid on

%% Wind plot
% figure(4)
% llll = plot(Rot_time,WindspeedX, Rot_time,WindspeedY);
% set(llll,'LineWidth',2);
% set(gca,'FontSize',18)
% %xlim([3e-4 1])
% %ylim([1e-10 1e-3])
% %set(gca,'YTick',10.^(-10:-3))
% xlabel('time (s)');
% ylabel('speed (mph)');
% %title('T240 seismometer');
% legend('X','Y',...
% 'Location','northeast')%'IPS X,RY noise')
% grid on

%% printing

if printFigs
figure(7)
FillPage('w')
%IDfig('K. Venkateswara')
saveas(gcf,'Earthquake Time Series 4-6-16.pdf')

figure(2)
FillPage('w')
% IDfig('K. Venkateswara')
saveas(gcf,'Earthquake ASD 4-6-16.pdf')

figure(11)
FillPage('w')
% IDfig('K. Venkateswara')
saveas(gcf,'Earthquake AngAvg 4-6-16.pdf')
figure(12)
FillPage('w')
% IDfig('K. Venkateswara')
saveas(gcf,'Earthquake VelAvg 4-6-16.pdf')
figure(13)
FillPage('w')
% IDfig('K. Venkateswara')
saveas(gcf,'Earthquake AngTim 4-6-16.pdf')
figure(14)
FillPage('w')
% IDfig('K. Venkateswara')
saveas(gcf,'Earthquake VelTim 4-6-16.pdf')
%
% figure(3)
% FillPage('w')
% IDfig('K. Venkateswara')
% saveas(gcf,'NoiseStudy2_Torque.pdf')
%
end