set(0,'DefaultAxesFontSize',16);
addpath(genpath('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/O1/H1/Scripts/CAL_PARAM/'));

% Channel names and analysis parameters; 
channelNames1 = {'H1:CAL-PCALY_TX_PD_OUT_DQ','H1:CAL-DELTAL_EXTERNAL_DQ'};
channelNames2 = {'H1:GDS-CALIB_STRAIN'};
segmentDuration = 1;
dataStartTime = 1237191970; % start time of the Pcal broadband injection done by Craig
dataDuration = 80;

% Get data
datastruct1 = get_data2(channelNames1, 'raw', dataStartTime, dataDuration);
datastruct2 = get_data2(channelNames2, 'raw', dataStartTime, dataDuration);

% Transfer function
windowFunct = hann(segmentDuration*datastruct1(1).rate);
noverlap = 0.5*segmentDuration*datastruct1(1).rate;
[deltal_pcal,freq] = tfestimate(datastruct1(1).data,datastruct1(2).data,windowFunct,noverlap,segmentDuration*datastruct1(1).rate,datastruct1(1).rate);
[gds_pcal,freq] = tfestimate(datastruct1(1).data,datastruct2.data,windowFunct,noverlap,segmentDuration*datastruct1(1).rate,datastruct1(1).rate);

% Dewhitening
fotonFileName = '/opt/rtcds/lho/h1/chans/H1CALCS.txt';
bankName = 'CS_DARM_DELTAL_EXTERNAL_WHITEN';
moduleNum = 2;
addpath('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/S7/Common/MatlabTools/');
white = readFilterResp(fotonFileName, bankName, moduleNum, freq);
rmpath('ligo/svncommon/CalSVN/aligocalibration/trunk/S7/Common/MatlabTools/');


deltal_pcal = deltal_pcal./white;
gds_pcal = gds_pcal*3995.4;

%%
figure(1)
subplot(2,1,1)
plot(freq,abs(deltal_pcal./gds_pcal),'b','LineWidth',2)
title('CAL DELTAL with whitening removed/ GDS STRAIN*3995.4')
xlim([20 200])
ylabel('Magnitude')
grid on
subplot(2,1,2)
plot(freq,rad2deg(angle(deltal_pcal./gds_pcal)),'b','LineWidth',2)
xlim([20 200])
xlabel('Frequency Hz')
ylabel('Phase (deg)')
grid on

print('CAL_DELTAL_over_GDS.png','-dpng')

figure(2)
subplot(2,1,1)
plot(freq,abs(deltal_pcal),'b','LineWidth',2)
title('CAL DELTAL with whitening removed/ PCAL Y TX PD')
xlim([20 200])
ylabel('Magnitude')
grid on
subplot(2,1,2)
plot(freq,rad2deg(angle(deltal_pcal)),'b','LineWidth',2)
xlim([20 200])
xlabel('Frequency Hz')
ylabel('Phase (deg)')
grid on
print('CAL_DELTAL_over_PCALY.png','-dpng')

figure(3)
subplot(2,1,1)
plot(freq,abs(gds_pcal),'b','LineWidth',2)
title('GDS/PCAL Y TX PD')
xlim([20 200])
ylabel('Magnitude')
grid on
subplot(2,1,2)
plot(freq,rad2deg(angle(deltal_pcal)),'b','LineWidth',2)
xlim([20 200])
xlabel('Frequency Hz')
ylabel('Phase (deg)')
grid on
print('GDS_over_PCALY.png','-dpng')

figure(4)
subplot(2,1,1)
plot(freq,abs(white),'b','LineWidth',2)
title('whitening filter')
xlim([20 200])
ylabel('Magnitude')
grid on
subplot(2,1,2)
plot(freq,rad2deg(angle(white)),'b','LineWidth',2)
xlim([20 200])
xlabel('Frequency Hz')
ylabel('Phase (deg)')
grid on

pcal_mean = mean(datastruct1(1).data);

figure(5)
plot(datastruct1(1).data)
title('PCALY TX PD')

