plothlts_matlabtfs.m

close all
clear all
clc

%% Hard-Coded Info
printFigs = 1;
saveData = 1;
printDiagDOFsOnly = 1;

ifo            = 'L1';           % ['H1' 'H2' 'L1' 'M1' 'X1' 'X2']
susID          = 'SR3';          % Optic Name or "QUAD" if on a test stand ['QUAD','ETMX','ITMY']
susNum         = '';             % Must be two digits, or empty string ['02','03','']
buildNum       = '';             % Must be two digits, or empty string ['01','02']
sagLevel       = 'M2';           % ['M0', 'R0', 'L1', 'L2', 'L3'] -- no damping anywhere else!
buildType      = 'hltsopt_metal';% ['wire','fiber']
meas.yymmdd    = '2016-06-15';   % e.g. 2012-04-23 for April 23rd 2012
meas.freqband  = '0p01to50Hz';
meas.time      = '1150064953';   % (whatever's after the date in the file name)
meas.sensCalib = true;           % true or false -- are sensors calibrated?
meas.author    = 'S. Aston';     % Let us know who did all the work!

isDamped       = 0;              % Plots the measured modeled closed loop against the measured one when set to 1.
% The filter file used is "dampingfilters_HSTS_20121023.mat"
% Name of pdf plot depends on this variable. If set to 1 damping_ON will be added to the filename, damping_OFF otherwise.

switch sagLevel
    case 'M1'
        dof.new = ['L';'T';'V';'R';'P';'Y'];
        dof.osem = ['T1';'T2';'T3';'LF';'RT';'SD'];
        nExcChans = 6;
        nRespChans = 12;
        nRespDofs = 6;

        % Matrix of Degrees of freedom which are expected to couple drive to
        % repsonse given the geometry of the HSTS. The matrix is intentionally
        % upper triangular, because it is assumed that transpose elements are
        % symmetric. To confirm this, the element and its trnaspose are plotted on
        % the same plot.

        expectedCoupling = [1  0  1  0  1  0;... % L
            0  1  0  1  0  0;... % T
            1  0  1  1  1  0;... % V
            0  1  1  1  1  0;... % R
            1  0  1  1  1  0;... % P
            0  0  0  0  0  1];    % Y

        meas.in.L = 1;
        meas.in.T = 2;
        meas.in.V = 3;
        meas.in.R = 4;
        meas.in.P = 5;
        meas.in.Y = 6;
        
        meas.out.osem.T1 = 1;
        meas.out.osem.T2 = 2;
        meas.out.osem.T3 = 3;
        meas.out.osem.LF = 4;
        meas.out.osem.RT = 5;
        meas.out.osem.SD = 6;
        meas.out.osem.L  = 7;
        meas.out.osem.T  = 8;
        meas.out.osem.V  = 9;
        meas.out.osem.R  = 10;
        meas.out.osem.P  = 11;
        meas.out.osem.Y  = 12;
        
    case {'M2','M3'}
        expectedCoupling = [1  1  1;... % L
                            1  1  1;... % P
                            1  1  1];   % Y
        
dof.new = ['L';'P';'Y'];
dof.osem = ['UL';'LL';'UR';'LR'];
nExcChans = 3;
nRespChans = 7;
nRespDofs = 3;

meas.in.L = 1;
meas.in.P = 2;
meas.in.Y = 3;
if sagLevel=='M2';
meas.out.osem.UL = 13;
meas.out.osem.LL = 14;
meas.out.osem.UR = 15;
meas.out.osem.LR = 16;
meas.out.osem.L  = 17;
meas.out.osem.P  = 18;
meas.out.osem.Y  = 19;
else
meas.out.osem.UL = 20;
meas.out.osem.LL = 21;
meas.out.osem.UR = 22;
meas.out.osem.LR = 23;
meas.out.osem.L  = 24;
meas.out.osem.P  = 25;
meas.out.osem.Y  = 26;    
end
end

svnDir = '/ligo/svncommon/SusSVN/sus/trunk/'; % LLO
% svnDir = '/ligo3/svncommon/SusSVN/sus/trunk/'; % LHO Test Stand Only

generalToolsDir = [svnDir 'Common/MatlabTools/'];
modelDir = [svnDir '/Common/MatlabTools/TripleModel_Production/'];

freqRange = [0.01 50]; % Hz
model.freq = logspace(log10(freqRange(1)), log10(freqRange(2)), 10000);


% Define DC scaling factor to calibrate data into
% (cts/ct) / (m/N) or (cts/ct) / (rad/N.m)
%% Calibrate data in rad/N.m or m/N

[calibration]=calib_hlts(susID,sagLevel,meas.sensCalib);

% Set Default Plot Parameters
set(0,'DefaultLineLineWidth',2) % Thick Lines
set(0,'DefaultAxesFontSize',16) % Nice big font
set(0,'DefaultAxesXGrid','on')  % Grid on by default
set(0,'DefaultAxesYGrid','on')
set(0,'DefaultAxesZGrid','on')

colorMatrix = [0.0 0.0 1.0;...
    0.0 0.5 1.0;...
    0.0 0.7 0.0;...
    0.0 1.0 0.0;...
    0.7 0.0 0.0;...
    1.0 0.0 0.0;...
    0.0 0.0 0.0;...
    0.5 0.5 0.5];

%% Generate Model

addpath(generalToolsDir);
addpath(modelDir);

[susModel,~,~,pendParams] = generate_Triple_Model_Production(model.freq,buildType,svnDir,true,isDamped);

%%
ls(svnDir)
if ~isempty(susNum)
    susTitle  = [ifo 'SUS' susID num2str(susNum) '-BUILD' buildNum];
    buildNum = num2str(buildNum);
    sagPath = [svnDir 'HLTS/' ifo '/' susID susNum '/BUILD' buildNum '/SAG' sagLevel '/'];
    dataDir = [sagPath 'Data/'];
    resultsDir = [sagPath 'Results/'];
else
    susTitle  = [ifo 'SUS' susID];
    sagPath = [svnDir 'HLTS/' ifo '/' susID '/SAG' sagLevel '/'];
    dataDir = [sagPath 'Data/'];
    resultsDir = [sagPath 'Results/'];
end

fileInfo = dir([dataDir meas.yymmdd '*' meas.time '*tf.mat']);
% fileInfo = dir([dataDir meas.yymmdd '*' meas.freqband '*.mat']);
% fileInfo = dir([dataDir meas.time '_*.mat']);
fileCell = struct2cell(fileInfo);  % added so 'strfind' can look through a cell
nFiles = length(fileInfo);

%%
disp('-----')
load([dataDir fileInfo.name]);
[data_full,meas.freq] = frdata(frd_all);

%% Expected Coupling
for iExcDOF = 1:nExcChans
    for iRespDOF = 1:nRespDofs
        if expectedCoupling(iExcDOF,iRespDOF)
            meas.eul2eul(iRespDOF,iExcDOF).f = squeeze(data_full(meas.out.osem.(dof.new(iRespDOF)),meas.in.(dof.new(iExcDOF)),:));
            if sagLevel=='M1'
                level='m1';
            elseif sagLevel=='M2'
                level='m2';
            else 
                level='m3';
            end
            if isDamped
                eval(['model.eul2eul(iRespDOF,iExcDOF).f = squeeze(susModel.dampedf(susModel.dampedout.' level '.disp.(dof.new(iRespDOF)),susModel.dampedin.' level '.drive.(dof.new(iExcDOF)),:));']);
            else
                eval(['model.eul2eul(iRespDOF,iExcDOF).f = squeeze(susModel.f(susModel.out.' level '.disp.(dof.new(iRespDOF)),susModel.in.' level '.drive.(dof.new(iExcDOF)),:));']);
            end
        else
            meas.eul2eul(iRespDOF,iExcDOF).f = zeros(length(meas.freq),1);
            model.eul2eul(iRespDOF,iExcDOF).f = zeros(length(model.freq),1);
        end
        meas.eul2eul(iRespDOF,iExcDOF).f = meas.eul2eul(iRespDOF,iExcDOF).f ./ (calibration);
    end
end

%% OSEM Basis
if sagLevel=='M1'
[OSEM2EUL, EUL2OSEM] = make_sushlts_projections;
invOSEM2EUL.(sagLevel) = inv(OSEM2EUL.(sagLevel));
invEUL2OSEM.(sagLevel) = inv(EUL2OSEM.(sagLevel));

% [oseminf,osemoutf,totalFilter,fotonFreq] = reconstructquadfotonfilters(meas.freq);

for iExcDOF = 1:nExcChans
    for iRespDOF = 1:nRespDofs
        meas.eul2osem(iRespDOF,iExcDOF).f = squeeze(data_full(meas.out.osem.(dof.osem(iRespDOF,:)),meas.in.(dof.new(iExcDOF)),:));
        meas.eul2osem(iRespDOF,iExcDOF).f = meas.eul2osem(iRespDOF,iExcDOF).f ./ (calibration);
        
        model.complete(iRespDOF,iExcDOF).f = squeeze(susModel.f(susModel.out.m1.disp.(dof.new(iRespDOF)),susModel.in.m1.drive.(dof.new(iExcDOF)),:));
    end
end

for iRow = 1:6
    for jColumn = 1:6
        model.eul2osem(iRow,jColumn).f = zeros(length(model.freq),1);
        model.osem2eul(iRow,jColumn).f = zeros(length(model.freq),1);
        model.osem2osem(iRow,jColumn).f = zeros(length(model.freq),1);
    end
end

for iRow = 1:6
    for jColumn = 1:6
        % Do the matrix multiplication manually, since the it doesn't work
        % for >2 D arrays (or for structures as I've done it)
        for kSummer = 1:nRespDofs
            model.eul2osem(iRow,jColumn).f = model.eul2osem(iRow,jColumn).f + invOSEM2EUL.(sagLevel)(iRow,kSummer) .* model.complete(kSummer,jColumn).f;
        end
        
        for kSummer = 1:nRespDofs
            model.osem2eul(iRow,jColumn).f = model.osem2eul(iRow,jColumn).f + model.complete(iRow,kSummer).f .* invEUL2OSEM.(sagLevel)(kSummer,jColumn);
        end
        for lSummer = 1:nRespDofs
            model.osem2osem(iRow,jColumn).f = model.osem2osem(iRow,jColumn).f + invOSEM2EUL.(sagLevel)(iRow,lSummer) .* model.osem2eul(lSummer,jColumn).f;
        end
    end
end
else
end
%%
iPlot = 1;
iPrintedPlot = 1;

for iExcDOF = 1:nExcChans
    for iRespDOF = 1:nRespDofs
        figure(iPlot)
        if iExcDOF == iRespDOF
            subplot(4,1,[1:3])
            ll=loglog(model.freq,abs(model.eul2eul(iRespDOF,iExcDOF).f),...
                meas.freq,abs(meas.eul2eul(iRespDOF,iExcDOF).f));
            title({susTitle;[sagLevel ' ' dof.new(iExcDOF) ' to ' dof.new(iRespDOF) ', ' meas.yymmdd]})
            %             set(ll(1),'Color',colorMatrix(1,:),'LineWidth',4)
            %             set(ll(2),'Color',colorMatrix(6,:))
            if strcmp(dof.new(iExcDOF),'L') || strcmp(dof.new(iExcDOF),'T') || strcmp(dof.new(iExcDOF),'V');
                ylabel('Magnitude (m/N)')
                ylim([1e-6 1])
                set(gca,'YTick',10.^(-6:0))
            elseif dof.new(iExcDOF)=='P'
                ylabel('Magnitude (rad/N.m)')
                ylim([1e-3 1e2])
                set(gca,'YTick',10.^(-3:2))
            else
                ylabel('Magnitude (rad/N.m)')
                ylim([1e-4 1e2])
                set(gca,'YTick',10.^(-4:2))
            end
            xlabel('Frequency (Hz)')
            xlim([1e-2 5e1])
            set(gca,'XTick',10.^(-2:2))
            set(gcf,'Position',[60 30 1050 750])
            text(0.05,0.9,0,['$\frac{Model}{Meas}$ @ 10mHz=' num2str(abs(model.eul2eul(iRespDOF,iExcDOF).f(1))/abs(meas.eul2eul(iRespDOF,iExcDOF).f(1)))],'interpreter','latex','Units','normalized','BackgroundColor',[.7 .9 .7],'FontSize',20)
            leg = legend(['Model (' buildType ')'],['Measurement / ' num2str(calibration) ],'Location','SouthWest');
            set(leg,'FontSize',16,'Interpreter','None')
            
            subplot(4,1,4)
            
            semilogx(model.freq,180/pi*angle(model.eul2eul(iRespDOF,iExcDOF).f),meas.freq,180/pi*angle(meas.eul2eul(iRespDOF,iExcDOF).f))
            xlabel('Frequency (Hz)')
            ylabel('Phase')
            xlim(freqRange)
            set(gca,'YTick',[-180:90:180])
            set(gca,'XTick',10.^(log10(freqRange(1)):log10(freqRange(2))))
            set(gcf,'Position',[60 30 1050 750])
            set(leg,'FontSize',16)
            
        elseif expectedCoupling(iExcDOF,iRespDOF) && iExcDOF > iRespDOF
            figure(iPlot)
            subplot(4,1,[1:3])
            ll=loglog(model.freq,abs(model.eul2eul(iRespDOF,iExcDOF).f),meas.freq,abs(meas.eul2eul(iRespDOF,iExcDOF).f),...
                model.freq,abs(model.eul2eul(iExcDOF,iRespDOF).f),meas.freq,abs(meas.eul2eul(iExcDOF,iRespDOF).f));
            title({susTitle;[sagLevel ' ' dof.new(iExcDOF) ' to ' dof.new(iRespDOF) ', ' meas.yymmdd]})
            set(ll([1 3]),'LineWidth',4);
%             for iColor = 1:8; set(ll(iColor),'Color',colorMatrix(iColor,:));end
            if strcmp(dof.new(iExcDOF),'L') || strcmp(dof.new(iExcDOF),'T') || strcmp(dof.new(iExcDOF),'V');
                ylabel('Magnitude (m/N, rad/N.m, or rad/N)')
                ylim([1e-6 1e2])
                set(gca,'YTick',10.^(-6:2))
            else
                ylabel('Magnitude (rad/N.m)')
                ylim([1e-5 1e2])
                set(gca,'YTick',10.^(-5:2))
            end
            xlabel('Frequency (Hz)')
            xlim(freqRange)
            set(gca,'XTick',10.^(log10(freqRange(1)):log10(freqRange(2))))
            set(gca,'XTick',10.^(-1:1))
            leg=legend(['Model (' dof.new(iExcDOF) ' to ' dof.new(iRespDOF) ')'],...
                ['Measurement / ' num2str(calibration) ' (' dof.new(iExcDOF) ' to ' dof.new(iRespDOF) ')'],...
                ['Model (' dof.new(iRespDOF) ' to ' dof.new(iExcDOF) ')'],...
                ['Measurement / ' num2str(calibration) ' (' dof.new(iRespDOF) ' to ' dof.new(iExcDOF) ')'],...
                'Location','SouthWest');
            set(leg,'FontSize',16)

            subplot(4,1,4)
            semilogx(model.freq,180/pi*angle(model.eul2eul(iRespDOF,iExcDOF).f),meas.freq,180/pi*angle(meas.eul2eul(iRespDOF,iExcDOF).f),...
                model.freq,180/pi*angle(model.eul2eul(iExcDOF,iRespDOF).f),meas.freq,180/pi*angle(meas.eul2eul(iExcDOF,iRespDOF).f))
            xlabel('Frequency (Hz)')
            ylabel('Phase')
            xlim(freqRange)
            set(gca,'YTick',[-180:90:180])
            set(gca,'XTick',10.^(log10(freqRange(1)):log10(freqRange(2))))
            set(gcf,'Position',[60 30 1050 750])
            set(leg,'FontSize',16)
        else
            close(iPlot)
        end
        % saveFileName = [resultsDir meas.yymmdd '_' susTitle '_' sagLevel]; % SMA - 17th Oct 2012 Added measurement time to filename
        saveFileName = [resultsDir meas.yymmdd '_' meas.time '_' susTitle '_' sagLevel];
        if printFigs && ((expectedCoupling(iExcDOF,iRespDOF) && iExcDOF > iRespDOF) || (iExcDOF == iRespDOF))
            plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(iExcDOF) '-' dof.new(iRespDOF) '_TF.pdf']};
            
            figure(iPlot)
            FillPage('w')
            IDfig(meas.author)
            saveas(gcf,plotName{iPrintedPlot,:})
            iPrintedPlot = iPrintedPlot + 1;
        end
        iPlot = iPlot + 1;
    end
end

if sagLevel=='M1'
%%
% EULER Output order : L T V R P Y = 1 2 3 4 5 6;
% OSEM Output order  : T1 T2 T3 LF RT SD = 1 2 3 4 5 6
% Input order        : L T V R P Y = 1 2 3 4 5 6;

figure(51)
ll = loglog(model.freq,abs(model.eul2eul(1,1).f),...
    meas.freq,abs(meas.eul2eul(1,1).f),...
    model.freq,abs(model.eul2eul(4,1).f),...
    meas.freq,abs(meas.eul2eul(4,1).f),...
    model.freq,abs(model.eul2osem(5,1).f),...
    meas.freq,abs(meas.eul2osem(5,1).f));
title({susTitle;[sagLevel ' ' dof.new(1) ' to ' dof.osem(4,:) dof.osem(5,:) ', ' meas.yymmdd]})
leg = legend('L (Model)',['L / ' num2str(calibration) ' (Meas)'],...
    'LF (Model)',['LF / ' num2str(calibration) ' (Meas)'],...
    'RT (Model)',['RT / ' num2str(calibration) ' (Meas)'],...
    'Location','SouthWest');
for iColor = 1:6; set(ll(iColor),'Color',colorMatrix(iColor,:));end
ylabel('Magnitude (m/N)')
xlabel('Frequency (Hz)')
xlim(freqRange)
ylim([1e-6 1e2])
set(gca,'YTick',10.^(-6:2),'XTick',10.^(-1:1))
set(ll([1 3 5]),'LineWidth',4)
set(leg,'FontSize',16);

%%
figure(52)
ll = loglog(model.freq,abs(model.eul2osem(6,2).f),...
    meas.freq,abs(meas.eul2osem(6,2).f));
title({susTitle;[sagLevel ' ' dof.new(2) ' to ' dof.osem(6,:) ', ' meas.yymmdd]})
leg = legend('SD (Model)',['SD / ' num2str(calibration) ' (Meas)'],...
    'Location','SouthWest');
for iColor = 1:2; set(ll(iColor),'Color',colorMatrix(iColor,:));end
ylabel('Magnitude (m/N)')
xlabel('Frequency (Hz)')
xlim(freqRange)
ylim([1e-6 1e2])
set(gca,'YTick',10.^(-6:2),'XTick',10.^(-1:1))
set(ll(1),'LineWidth',4)
set(leg,'FontSize',16);

%%
figure(53)
ll = loglog(model.freq,abs(model.eul2eul(3,3).f),...
    meas.freq,abs(meas.eul2eul(3,3).f),...
    model.freq,abs(model.eul2osem(1,3).f),...
    meas.freq,abs(meas.eul2osem(1,3).f),...
    model.freq,abs(model.eul2osem(2,3).f),...
    meas.freq,abs(meas.eul2osem(2,3).f),...
    model.freq,abs(model.eul2osem(3,3).f),...
    meas.freq,abs(meas.eul2osem(3,3).f));
title({susTitle;[sagLevel ' ' dof.new(3) ' to ' dof.osem(1,:) dof.osem(2,:) dof.osem(3,:) ', ' meas.yymmdd]})
leg = legend('V (Model)',['V / ' num2str(calibration) ' (Meas)'],...
    'T1 (Model)',['T1 / ' num2str(calibration) ' (Meas)'],...
    'T2 (Model)',['T2 / ' num2str(calibration) ' (Meas)'],...
    'T3 (Model)',['T3 / ' num2str(calibration) ' (Meas)'],...
    'Location','SouthWest');
for iColor = 1:8; set(ll(iColor),'Color',colorMatrix(iColor,:));end
ylabel('Magnitude (m/N)')
xlabel('Frequency (Hz)')
xlim(freqRange)
ylim([1e-6 1e2])
set(gca,'YTick',10.^(-6:2),'XTick',10.^(-1:1))
set(ll([1 3 5 7]),'LineWidth',4)
set(leg,'FontSize',16);

%%
figure(54)
ll = loglog(model.freq,abs(model.eul2eul(4,4).f),...
    meas.freq,abs(meas.eul2eul(4,4).f),...
    model.freq,abs(model.eul2osem(1,4).f),...
    meas.freq,abs(meas.eul2osem(1,4).f),...
    model.freq,abs(model.eul2osem(2,4).f),...
    meas.freq,abs(meas.eul2osem(2,4).f),...
    model.freq,abs(model.eul2osem(3,4).f),...
    meas.freq,abs(meas.eul2osem(3,4).f));
title({susTitle;[sagLevel ' ' dof.new(4) ' to ' dof.osem(1,:) dof.osem(2,:) dof.osem(3,:) ', ' meas.yymmdd]})
leg = legend('R (Model)',['R / ' num2str(calibration) ' (Meas)'],...
    'T1 (Model)',['T1 / ' num2str(calibration) ' (Meas)'],...
    'T2 (Model)',['T2 / ' num2str(calibration) ' (Meas)'],...
    'T3 (Model)',['T3 / ' num2str(calibration) ' (Meas)'],...
    'Location','SouthWest');
for iColor = 1:8; set(ll(iColor),'Color',colorMatrix(iColor,:));end
ylabel('Magnitude (rad/N.m or m/N.m)')
xlabel('Frequency (Hz)')
xlim(freqRange)
ylim([1e-6 1e2])
set(gca,'YTick',10.^(-6:2),'XTick',10.^(-1:1))
set(ll([1 3 5 7]),'LineWidth',4)
set(leg,'FontSize',16);

%%
figure(55)
ll = loglog(model.freq,abs(model.eul2eul(5,5).f),...
    meas.freq,abs(meas.eul2eul(5,5).f),...
    model.freq,abs(model.eul2osem(2,5).f),...
    meas.freq,abs(meas.eul2osem(2,5).f),...
    model.freq,abs(model.eul2osem(3,5).f),...
    meas.freq,abs(meas.eul2osem(3,5).f));
title({susTitle;[sagLevel ' ' dof.new(5) ' to ' dof.osem(2,:) dof.osem(3,:) ', ' meas.yymmdd]})
leg = legend('P (Model)',['P / ' num2str(calibration) ' (Meas)'],...
    'T2 (Model)',['T2 / ' num2str(calibration) ' (Meas)'],...
    'T3 (Model)',['T3 / ' num2str(calibration) ' (Meas)'],...
    'Location','SouthWest');
for iColor = 1:6; set(ll(iColor),'Color',colorMatrix(iColor,:));end
ylabel('Magnitude (rad/N.m or m/N.m)')
xlabel('Frequency (Hz)')
xlim(freqRange)
ylim([1e-6 1e2])
set(gca,'YTick',10.^(-6:2),'XTick',10.^(-1:1))
set(ll([1 3 5 ]),'LineWidth',4)
set(leg,'FontSize',16);

%%
figure(56)
ll = loglog(model.freq,abs(model.eul2eul(6,6).f),...
    meas.freq,abs(meas.eul2eul(6,6).f),...
    model.freq,abs(model.eul2osem(4,6).f),...
    meas.freq,abs(meas.eul2osem(4,6).f),...
    model.freq,abs(model.eul2osem(5,6).f),...
    meas.freq,abs(meas.eul2osem(5,6).f));
title({susTitle;[sagLevel ' ' dof.new(6) ' to ' dof.osem(4,:) dof.osem(5,:) ', ' meas.yymmdd]})
leg = legend('Y (Model)',['Y / ' num2str(calibration) ' (Meas)'],...
    'LF (Model)',['LF / ' num2str(calibration) ' (Meas)'],...
    'RT (Model)',['RT / ' num2str(calibration) ' (Meas)'],...
    'Location','SouthWest');
for iColor = 1:6; set(ll(iColor),'Color',colorMatrix(iColor,:));end
ylabel('Magnitude (rad/N.m or m/N.m)')
xlabel('Frequency (Hz)')
xlim(freqRange)
ylim([1e-6 1e2])
set(gca,'YTick',10.^(-6:2),'XTick',10.^(-1:1))
set(ll([1 3 5]),'LineWidth',4)
set(leg,'FontSize',16);

%%
if printFigs
    figure(51)
    FillPage('w')
    IDfig(meas.author)
    plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(1) '-' dof.osem(2,:) dof.osem(3,:) '_TF.pdf']};
    saveas(gcf,plotName{iPrintedPlot,:})
    iPrintedPlot = iPrintedPlot + 1;
    
    figure(52)
    FillPage('w')
    IDfig(meas.author)
    plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(2) '-' dof.osem(6,:) '_TF.pdf']};
    saveas(gcf,plotName{iPrintedPlot,:})
    iPrintedPlot = iPrintedPlot + 1;
    
    figure(53)
    FillPage('w')
    IDfig(meas.author)
    plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(3) '-' dof.osem(4,:) dof.osem(5,:) '_TF.pdf']};
    saveas(gcf,plotName{iPrintedPlot,:})
    iPrintedPlot = iPrintedPlot + 1;
    
    figure(54)
    FillPage('w')
    IDfig(meas.author)
    plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(4) '-' dof.osem(4,:) dof.osem(5,:) '_TF.pdf']};
    saveas(gcf,plotName{iPrintedPlot,:})
    iPrintedPlot = iPrintedPlot + 1;
    
    figure(55)
    FillPage('w')
    IDfig(meas.author)
    plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(5) '-' dof.osem(1,:) dof.osem(2,:) dof.osem(3,:) '_TF.pdf']};
    saveas(gcf,plotName{iPrintedPlot,:})
    iPrintedPlot = iPrintedPlot + 1;
    
    figure(56)
    FillPage('w')
    IDfig(meas.author)
    plotName(iPrintedPlot,:) = {[saveFileName '_' dof.new(6) '-' dof.osem(2,:) dof.osem(3,:) '_TF.pdf']};
    saveas(gcf,plotName{iPrintedPlot,:})
    iPrintedPlot = iPrintedPlot + 1;
end
else
end
%% Make Summary .pdf of all figures
if isDamped==true;
    meas.damping='ON';
else
    meas.damping='OFF';
end

% [diag TFs]       [off-diag TFs]   [osembasis TFs]
mergeOrder = [1 2 4 7 11 12    3 5 6 8 9 10     13 14 15 16 17 18];
mergeOrderDOF = [1 2 4 7 11 12]; % Diag TF only
if strcmp(sagLevel,'M2') || strcmp(sagLevel,'M3')
    mergeOrder = [1 3 6 2 4 5];
    mergeOrderDOF = [1 3 6];
end
if printFigs
    mergeCommand = ['pdfmerge ' saveFileName '_damp_' meas.damping '_ALL_TFs.pdf'];
    if printDiagDOFsOnly == 1;
        for iFile = mergeOrderDOF
            mergeCommand = [mergeCommand ' ' plotName{iFile}];
        end
    else
        for iFile = mergeOrder
            mergeCommand = [mergeCommand ' ' plotName{iFile}];
        end
    end
        
    % SMA - 01/27/2017 - Added pdfmerge workaround for SL7 (https://services.ligo-la.caltech.edu/FRS/show_bug.cgi?id=7240)
    savedPath = getenv('LD_LIBRARY_PATH');
    setenv('LD_LIBRARY_PATH', '');
    system(mergeCommand);
    setenv('LD_LIBRARY_PATH', savedPath);

    disp(['saved ' saveFileName '_damp_' meas.damping '_ALL_TFs.pdf'])
    
    for iFile = mergeOrder
        system(['rm ' plotName{iFile}]);
    end
end

%% Save the Data

if saveData
    save([saveFileName '.mat'],'meas','model')
end