
% This tagged model uses damping filters from 2012-10-23, which 
% were a replica of L1 SUS MC2's damping filters at the time. 
% See LLO aLOG 6041 for comparison of TOP to TOP TFs model vs. measurement
% to get a feel for the accuracy.
% https://alog.ligo-la.caltech.edu/aLOG/index.php?callRep=6041

% The tagged model can be found here:
% https://redoubt.ligo-wa.caltech.edu/svn/sus/trunk/Common/SusModelTags/Matlab/
% hstsmodelproduction-rev4067_ssmake3MBf-rev1891_hstsopt_metal-rev2039_released-2013-01-31.mat

%% Change all this stuff to point to where you've downloaded the model
modelDir = '/ligo/svncommon/SusSVN/sus/trunk/Common/SusModelTags/Matlab/';
modelTag = 'hstsmodelproduction-rev4067_ssmake3MBf-rev1891_hstsopt_metal-rev2039_released-2013-01-31.mat';

dataDir = '/ligo/svncommon/SusSVN/sus/trunk/HSTS/H1/PRM/SAGM3/Data/';
fileNames = {'2018-02-13_2350_H1SUSPRM_M3_WhiteNoise_L_0p01to50Hz_tf.txt';...
             '2018-02-13_2350_H1SUSPRM_M3_WhiteNoise_P_0p01to50Hz_tf.txt';...
             '2018-02-13_2350_H1SUSPRM_M3_WhiteNoise_Y_0p01to50Hz_tf.txt'};

%% Change this to false if you don't need to print the figure
printFigs = true;
resultsDir = '/Users/kissel/Desktop/scratch/';
figTag = '2018-02-13_SUSHSTS_M3toM3_LandP_Model_vs_Meas_Responses.pdf';

%% Load in the Model
load([modelDir modelTag]);

%% Calibrate the Data
sagLevel = 3;
dac_drive=20/2^18; %V/cts
coil_driver=[0.011919 0.0028 0.003]; %A/V G1100968-v6
coil_magnet=[0.963 0.0158 0.00281]; %N/A
osem_sensitivity=0.7/76.29; %mm/uA
sat_amp=1/240e3; %A/V
adc=40/2^16; %V/cts
calibration = (0.4)*1e6*dac_drive*coil_driver(sagLevel)*coil_magnet(sagLevel);

%% Load the Data
excL  = 1; excP  = 2; excY  = 3;
respL = 5; respP = 6; respY = 7;

for iExcDOF = 1:3
    data(iExcDOF).raw = load([dataDir fileNames{iExcDOF}]);
    data(iExcDOF).freq = data(iExcDOF).raw(:,1);
    
    for iRespChan = 1:7
        realInd = 2*iRespChan;
        imagInd = 2*iRespChan + 1;
        
        data(iExcDOF).tf(:,iRespChan) = (data(iExcDOF).raw(:,realInd) + sqrt(-1)*data(iExcDOF).raw(:,imagInd))./calibration;
    end
end



%% Plot the desired transfer functions
figure(1);
   subplot(221)
   loglog(freq,abs(squeeze(hstsModel.dampedf(hstsModel.dampedout.m3.disp.L,hstsModel.dampedin.m3.drive.L,:))),...
          data(excL).freq,abs(data(excL).tf(:,respL)));
   xlabel('Frequency / Hz'); xlim([1e-1 1e1]); set(gca,'XTick',10.^(-1:1))
   ylabel('Magnitude / m.N^{-1}'); ylim([1e-6 1e0]); set(gca,'YTick',10.^(-6:0))
   title('M3 L Drive to M3 L Response')
   legend('Model','2018-02-13 H1SUSPRM Meas','Location','SouthWest')
   
   subplot(223)
   loglog(freq,abs(squeeze(hstsModel.dampedf(hstsModel.dampedout.m3.disp.L,hstsModel.dampedin.m3.drive.P,:))),...
          data(excL).freq,abs(data(excL).tf(:,respP)));
   xlabel('Frequency / Hz'); xlim([1e-1 1e1]); set(gca,'XTick',10.^(-1:1))
   ylabel('Magnitude / m.(N.m)^{-1}'); ylim([1e-6 2e0]); set(gca,'YTick',10.^(-6:0))
   title('M3 P Drive to M3 L Response')
   
   subplot(222)
   loglog(freq,abs(squeeze(hstsModel.dampedf(hstsModel.dampedout.m3.disp.P,hstsModel.dampedin.m3.drive.L,:))),...
          data(excP).freq,abs(data(excP).tf(:,respL)));
   xlabel('Frequency / Hz'); xlim([1e-1 1e1]); set(gca,'XTick',10.^(-1:1))
   ylabel('Magnitude / rad.N^{-1}'); ylim([1e-6 1e1]); set(gca,'YTick',10.^(-5:1))
   title('M3 L Drive to M3 P Response')
   
   subplot(224)
   loglog(freq,abs(squeeze(hstsModel.dampedf(hstsModel.dampedout.m3.disp.P,hstsModel.dampedin.m3.drive.P,:))),...
          data(excP).freq,abs(data(excP).tf(:,respP)));
   xlabel('Frequency / Hz'); xlim([1e-1 1e1]); set(gca,'XTick',10.^(-1:1))
   ylabel('Magnitude / rad.(N.m)^{-1}'); ylim([1e-3 2e2]); set(gca,'YTick',10.^(-3:2))
   title('M3 P Drive to M3 P Response')
  
   
%% Print the figure, if need be   
if printFigs
    figure(1)
    FillPage('w')
    IDfig
    saveas(gcf,[resultsDir figTag])
end