clear all

%DC_photocurrent=20; %mA
%optical_gain=3.3;%ma/pm

%DC_offset= 1e-12*2*DC_photocurrent/optical_gain; % 12 picometers

%fund_curr=0.54;%mA/rtHz
%harm_curr=4.61e-4; %mA/rt Hz 12 Hz line
%BW=0.03; %BW used in DTT

%fund_amp=fund_curr*sqrt(BW);
%harm_amp=harm_curr*sqrt(BW);

%ratio=harm_amp/fund_amp

%predicted_ratio=fund_amp/(4*DC_photocurrent)

%check against measured current

addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/
% 

start =[1137307997,1137631136,1137630496]; %6:53 UTC Jan 20 2015
dt = 2*60;
gps = start;
chan = {'H1:OMC-DCPD_SUM_OUT_DQ','H1:LSC-ASAIR_A_RF45_Q_ERR_DQ'};%, 'H1:LSC-ASAIR_A_RF45_I_ERR_DQ'};'H1:LSC-DARM_IN1_DQ',
gwd = GWData();
gwd.make_kerberos_ready;
gwd.site_info(2).server = 'nds2.ligo-wa.caltech.edu';
gwd.site_info(2).port = 31200;
for ii=1:length(start)
    [data(ii,:,:), t, info] = gwd.fetch(start(ii), dt, chan);
end

%save('AS45_vs_DCPDs_DARM_boost_data.mat', 'data')

data_scaled=data;
AS45_scaled=squeeze(data_scaled(:,:,2)./(6.81e-7/0.52));
AS45_means=mean(AS45_scaled,2);
for ii=1:length(start)
   AS45_detrend(ii,:)=AS45_scaled(ii,:)-AS45_means(ii); 
end

%make a filter to apply
Fs=16384;  %sampling frequency
Wn=5/(Fs/2);  %you cut off frequency, normalized to the myquist frequency
[z,p,k]=butter(6,Wn);
for ii=1:length(start)
    AS45_filtered(ii,:)=apply_filter(z,p,k,AS45_detrend(ii,:), 100);
end

%find the sum of the DC PD residual with a lock point offset residual. 
DCPDs=squeeze(data_scaled(:,:,1));
DC_term=mean(DCPDs,2);
for ii=1:length(start)
    res(:,ii)=DCPDs(ii,:)-DC_term(ii);
    res_w_err(:,ii)=DCPDs(ii,:)-DC_term(ii)+AS45_filtered(ii,:);
    upconvo(:,ii)=(res(:,ii)).^2/(4*DC_term(ii));
    upconvo_w_err(:,ii)=(res_w_err(:,ii)).^2/(4*DC_term(ii));
end
%make spectra
T=1/Fs;    %sample time
stepSize = 30; %1 second
overlap = 0.5;
NFFT = stepSize*Fs;
w = hanning(NFFT);
for ii=1:length(start)
    %scale AS45, there is an 80 in the input matrix when we are locked on
    %RF, power norm matrix has TRX (1254 cnts) and ain of 0.37
    %fudge factor of 3.6 is just chose to match peak at 2.8 Hz
    data_scaled(ii,:,2)=data(ii,:,2)/3.97;%*80/((sqrt(1254)*0.37));
    [P_DCPD(:,ii), f]=pwelch(data_scaled(ii,:,1),w,overlap*NFFT,NFFT,Fs);
    [P_AS45(:,ii), f]=pwelch(AS45_detrend(ii,:),w,overlap*NFFT,NFFT,Fs);
    [P_AS45filt(:,ii), f]=pwelch(AS45_filtered(ii,:),w,overlap*NFFT,NFFT,Fs);
    [P_upconvo(:,ii), f]=pwelch(upconvo(:,ii),w,overlap*NFFT,NFFT,Fs);
    [P_upconvo_w_err(:,ii), f]=pwelch(upconvo_w_err(:,ii),w,overlap*NFFT,NFFT,Fs);
    [P_res(:,ii), f]=pwelch(res(:,ii),w,overlap*NFFT,NFFT,Fs);
    [P_res_w_err(:,ii), f]=pwelch(res_w_err(:,ii),w,overlap*NFFT,NFFT,Fs);
end

h1=figure(1);
loglog(f', sqrt(P_DCPD(:,1)),'b',f',sqrt(P_AS45(:,1)),'--b',...
    f', sqrt(P_DCPD(:,2)),'r',f',sqrt(P_AS45(:,2)),'--r',...
    f', sqrt(P_DCPD(:,3)),'k',f',sqrt(P_AS45(:,3)),'--k',...
    'LineWidth',2)
xlim([-Inf,30])
%ylim([1e-10,1])
ylabel('mA rt Hz on DCPDs')
xlabel('Frequency')
grid on
legend('DCPD photocurrent Jan 20 with 6 Hz excitation','AS45 Jan 20',...
    'DCPD photocurrent after feedforward tunning', 'AS45 after feedforward tunning',...
    'DCPD photocurrent FF tuning and DARM boost', 'AS45 after FF and DARM boost')
%print(h1, '-dpng','DC_readout_upconvo_ff_tuning.png')

h2=figure(2);
loglog(f', sqrt(P_DCPD(:,2)),'r',f',sqrt(P_AS45(:,2)),'--r',...
    f', sqrt(P_DCPD(:,3)),'k',f',sqrt(P_AS45(:,3)),'--k',...
        f', sqrt(P_AS45filt(:,2)),'*r',f',sqrt(P_AS45filt(:,3)),'*k',...
    'LineWidth',2)
xlim([-Inf,1000])
%ylim([1e-10,1])
ylabel('mA rt Hz on DCPDs')
xlabel('Frequency')
grid on
legend('DCPD photocurrent after feedforward tunning', 'AS45 after feedforward tunning',...
    'DCPD photocurrent FF tuning and DARM boost', 'AS45 after FF and DARM boost', ...
    'Filtered AS45 without boost', 'Filtered AS45 with boost')

h4=figure(4);
loglog(f', sqrt(P_DCPD(:,2)),'r',f',sqrt(P_AS45(:,2)),'--r',...
    f', sqrt(P_DCPD(:,3)),'k',f',sqrt(P_AS45(:,3)),'--k',...
    'LineWidth',2)
xlim([-Inf,1000])
%ylim([1e-10,1])
ylabel('mA rt Hz on DCPDs')
xlabel('Frequency')
grid on
legend('DCPD photocurrent nominal DARM loop', 'AS45 nominal DARM loop',...
    'DCPD photocurrent DARM boost and resG', 'AS45 after DARM boost and resG')
print(h4, 'AS45_vs_DCPDs_DARM_boost.png')


h3=figure(3);
loglog(f', sqrt(P_res(:,2)),'r',f', sqrt(P_res_w_err(:,2)), 'b',...
    f', sqrt(P_upconvo(:,2)), 'k',f', sqrt(P_upconvo_w_err(:,2)),'m',...
    'LineWidth',2)
xlim([-Inf,1000])
%ylim([1e-10,1])
ylabel('mA rt Hz on DCPDs')
xlabel('Frequency')
grid on
legend('DCPD residual', 'residual with lock point errors seen by AS45 added', ...
    'upconversion DPCD residual only', 'upconversion if AS45 fluctuations are DARM offset fluctuations')
print(h3, 'AS45_lock_point_errors_upconvo.png')