addpath('/ligo/svncommon/CalSVN/aligocalibration/trunk/Common/MatlabTools/utils');
addpath('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/O2/DARMmodel/src');

% load the parameter files
IFOindepPars = '/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/ER10/Common/params/IFOindepParams.conf';

IFOdepPars   = {'/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/ER10/H1/params/H1params.conf'};

IFOmeasPars  = {'/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/ER10/H1/params/2016-11-12/H1params_2016-11-12.conf'};

extract = {'calcsResid2gdsInvSensOut','pcalCorrection'};

f = [1:0.1:16384]';
overwrite = [];
[~,~,~,~,DARMvars] = computeDARM(IFOindepPars, IFOdepPars{1}, IFOmeasPars{1}, f, overwrite, extract);

dataCal = load('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/PreER10/H1/Scripts/ControlRoomCalib/pcal2darm_calib.txt');

sys = zpk(-2*pi*[30 30 30 30 30 30], -2*pi*[0.3 0.3 0.3 0.3 0.3 0.3], 1e-12);
dewhitenEXT = squeeze(freqresp(sys, 2*pi*f));

EXT_Pcal_correct = -dewhitenEXT./DARMvars{1,2}./DARMvars{2,2}.*exp(sqrt(-1)*2*pi*f*61e-6);

dataCal1 = interp1(dataCal(:,1), db2mag(dataCal(:,2)).*exp(sqrt(-1)*dataCal(:,3)*2*pi/360), f);

figure
subplot(2,2,1)
loglog(f, abs(dataCal1))
hold on
loglog(f, abs(EXT_Pcal_correct))
grid on
legend('old','new')
subplot(2,2,3)
semilogx(f, rad2deg(angle(dataCal1)))
hold on
semilogx(f, rad2deg(angle(EXT_Pcal_correct)))
xlabel('Freq (Hz)')
grid on
legend('old','new')
subplot(2,2,2)
semilogx(f, abs(EXT_Pcal_correct./dataCal1))
grid on
ylabel('abs[New/Old]')
subplot(2,2,4)
semilogx(f, rad2deg(angle(EXT_Pcal_correct./dataCal1)))
grid on
ylabel('Phase difference [New/Old]')
xlabel('Freq (Hz)')
