% find time when we were adjusting modulation depths to check ratio of
% powers in carrier and sidebands at different ports. 

%% get the data
addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/

slow_channels = {'H1:LSC-MOD_RF9_AM_RFOUT','H1:LSC-MOD_RF45_AM_RFOUT',...
        'H1:ASC-OMC_NSUM_OUTMON',...
    'H1:ASC-POP_A_NSUM_OUT16', 'H1:ASC-POP_B_NSUM_OUT16',... %    'H1:LSC-POP_A_LP_OUT16', LSC POP has more noise
    'H1:LSC-REFL_A_LF_OUT16', 'H1:LSC-REFL_B_LF_OUT16',...
    'H1:ASC-X_PWR_CIRC_OUT16','H1:ASC-Y_PWR_CIRC_OUT16',...
};


% fetch channel data
gwd = GWData();
gwd.make_kerberos_ready;
gwd.site_info(2).server = 'nds.ligo.caltech.edu';
gwd.site_info(2).port = 31200;
start_time = 1239497188;
dt = 300;
[slow_data, t, ~] = gwd.fetch(start_time, dt, slow_channels);


%we want to look at the sum of POP A+B, sum of refl A+B, 
data = [slow_data(:,1:3), slow_data(:,4)+slow_data(:,5), slow_data(:,6)+slow_data(:,7), slow_data(:,8)+slow_data(:,9)];
%only do this once each time you fetch data
%note the indeces so you don't have to keep track
dB9 = 1;
dB45 = 2;
as = 3;
pop = 4;
refl = 5;
tr = 6;
%%

means = mean(data(1:100, :)); %normalize to the first 100 points, these are the first few seconds (nominal powers)
for ii = 1:size(data,2)
    norm_data(:,ii) = data(:,ii)./means(ii);
end

figure(1)
subplot(3,1,1)
plot(t, slow_data(:,1),t, slow_data(:,2))
legend('9MHz level', 'RF45 level')
ylabel('dBm')
subplot(3,1,2)
plot(t, slow_data(:,4),t, slow_data(:,5))%
legend('ASC pop A', 'ASC pop B', 'A+B')
legend('POP A LF')
subplot(3,1,3)
plot(t, slow_data(:,6),t, slow_data(:,7))
legend('Refl A LF', 'REFL B LF',  'REFL A+B')

figure(2)
subplot(2,1,1)
plot(t, slow_data(:,1),t, slow_data(:,2))
legend('9MHz level', 'RF45 level')
ylabel('dBm')
subplot(2,1,2)
plot(t, data(:,3),t,data(:,4), t, data(:,5), t, data(:,6))
legend('AS', 'POP', 'REFL', 'TR')

%find times at nominal, max 45, and max 9
start_nom = 0;
end_nom = 42; %we were at nominal settings for first 42 seconds of this time
start_45 = 65.5;
end_45 = 88.3;
start_9 = 194.3;
end_9 = 210;

time_vector = [start_nom, end_nom, start_45, end_45, start_9, end_9];
for ii=1:length(time_vector)
    time_index(ii) = find(t>time_vector(ii),1);
end
%take means of the data for each of the periods
means_nom = mean(data(time_index(1):time_index(2),:));
means_45 = mean(data(time_index(3):time_index(4),:));
means_9 = mean(data(time_index(5):time_index(6),:));

%gamma ratio
mod_change_9 =  10^((means_9(dB9)-means_nom(dB9))/20);
mod_change_45 =  10^((means_45(dB45)-means_nom(dB45))/20);

%you can estimate modulation depth using the drop in arm power when the mod
%depth is changed (alignment changes could mess it up).
tr_drop_9 = means_9(tr)./means_nom(tr);
tr_drop_45 = means_45(tr)./means_nom(tr);

% gamma45 =0.1352;
% (besselj(0,gamma45*mod_change_45)/besselj(0,gamma45))^2
car_drop = @(gamma, mod_change)(besselj(0,gamma*mod_change)/besselj(0,gamma))^2;
solve_for_gamma45 = @(gamma45)car_drop(gamma45,mod_change_45)-tr_drop_45;
gamma45_infered = fzero(solve_for_gamma45,0.172);  %inital gues of 0.172 for gamma 45 (alog 48513)

solve_for_gamma9 = @(gamma9)car_drop(gamma9,mod_change_9)-tr_drop_9;
gamma9_infered = fzero(solve_for_gamma9,0.135);  %inital gues
%you can see from the normalized POP build ups that when we change the 9
%MHz modulation depth, we change some alignments

%instead of using these infered gammas, use our best gueses 
gamma9 = 0.135;
gamma45 = 0.172;

%now find the amount of power in each field at each detector
%make some vectors of powers at each port for different measurements
m = [means_nom(as:tr);...
    means_45(as:tr);...
    means_9(as:tr)]; %each column is a different port, first row is nominal modulation depths, 2nd row is 45 increased, 3rd row is 9 increased

%here is a matrix A, the equation A[Pc;P45;P9] = m(column) for each column
%of m gives Pc/Jo^2Jo^2 is power in the carrier at that diode, P45/J1gamma45^2 + P9/J1gamma9^2 are the
%sideband powers all for the nominal modulation depth
A = [besselj(0,gamma45)^2*besselj(0,gamma9)^2, besselj(1,gamma45)^2, besselj(1,gamma9)^2;...
    besselj(0,mod_change_45*gamma45)^2*besselj(0,gamma9)^2, besselj(1,mod_change_45*gamma45)^2, besselj(1,gamma9)^2;...
    besselj(0,gamma45)^2*besselj(0,mod_change_9*gamma9)^2, besselj(1,gamma45)^2, besselj(1,mod_change_9*gamma9)^2];

for ii = 1:4
    %1 is as, 2 is pop, 3 is refl, 4 is tr
    p(:,ii)= A^-1*m(:,ii);
end

%P needs to be multipied byt the besselJ factors to recover power at
%nominal modulation depths
%first row is for carrier
P(1,:) = p(1,:)*besselj(0,gamma45)^2*besselj(0,gamma9)^2;
P(2,:) = p(2,:)*besselj(1,gamma45)^2;
P(3,:) = p(3,:)*besselj(1,gamma9)^2;

tot=sum(P); %this is the total power on each diode, should be the same as means_nom(as:tr)

%find the ratio of different powers at different ports
for ii=1:4
    percent(:,ii)=100*P(:,ii)/tot(ii);
end
